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ABSTRACT 


The statistical properties of low-1eve1 wind-turbulence data 
were obtained with the model 1080 total vector anemometer and the 
model 1296 dual split-film anemometer, both manufactured by Thermo- 
Systems Incorporated. The data obtained from the above fast-response 
probes were compared with the results obtained from a pair of Gill 
propeller anemometers. 

The digitized time series representing the three velocity com- 
ponents and the temperature were each divided into a number of blocks, 
the length of which depended on the lowest frequency of Interest and 
also on the storage capacity of the available computer. A moving- 
average and differencing high-pass filter was used to remove the trend 
and the low frequency components in the time series. 

The degree of nonstationarity of each time series was determined 
by using a non-parametric statistical test on the statistical quantities 
calculated for each block in the time series. Besides the mean, the. 
variances and the covariances of the fluctuating velocity components 
and the fluctuating temperature, spectral and cross-spectral estimates of 
each of the time series were obtained with the use of the fast Fourier 
transformation (F.F.T) technique. A special F.F.T. algorithmic with a 
no-bit reversal procedure for the analysis of series of long duration 
and high sampling rate (200 samples per second) has been developed. 

In addition, a time series representing the streamwise fluctuating 
velocity component was simulated from the semi-empirical von Karman spectrum 
equation. The spectrum calculated from the simulated time series was com- 
pared with the original spectrum function from which the data was obtained. 

The calculated results for each of the anemometers used are represented 
in graphical or tabulated form. The Fortran program for the entire data 
analysis procedure is listed in appendix B. 
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CHAPTER I 
INTRODUCTION 

With the increasing necessity of understanding the effect of the 
turbulence on human activities in the atmospheric boundary layer, the 
knowledge of wind and temperature fluctuations becomes essential. 

This dissertation presents a statistical analysis of fluctuating wind 
data consisting of three velocity components for each of the three 
spatial directions and temperature in the time domain. Correlations 
and spectra will be calculated by using the experimental data ob- 
tained from T.S.I. split - film anemometers and a consequent data 
acquisition system as reported by Tieleman et. al. [95]. 

1.1 STATEMENT OF THE PROBLEM 

Although geophysical flows provide much larger Reynolds numbers 
than wind tunnel flows, relatively few accurate results of the wind 
turbulence measurements in the atmospheric boundary layer are avail- 
able because of the complexity and inadequacy of the required instru- 
mentation. Existing measuring equipment usually cannot provide 
enough information for the three dimensional structure of the turbu- 
lence and its evolution with time. The selection of the three dimen- 
sional split-film anemometer (TSI-1080D) in this research program has 
many advantages over most of the presently used anemometers [95]. 
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Insufficient knowledge of the ground wind structure is due to 
a lack of extensive experiments and due to large volumes of data 
necessary to obtain accurate and meaningful results. In order to 
understand the microstructure of turbulence, we must not only consider 
the energy spectrum but also several statistical quantities which are 
relevant to the basic mechanics of the motion of the air. The measure 
ment of the ground wind structure and its statistical analysis becomes 
essential to the solution of such practical problems in the atmo- 
spheric boundary layer as low altitude operations of aircraft, design 
of many engineering structures, and atmospheric diffusion. It is 
therefore imperative to improve statistical calculation techniques of 
the measured quantities especially when dealing with large volumes of 
data even beyond the storage capacity of the existing computers. 

The purpose of this dissertation is to analyse the data based up- 
on the general aspects of time series analysis. First, the analysis 
in the time domain where the means, variances and covariances of 
blocked subseries are calculated. Second, the analysis in the fre- 
quency domain where the power- and cross- spectral density estimations 
are calculated. 

The power spectrum which provides an insight to the spectral dis- 
tribution of the kinetic energy of the turbulence, is basic to the 
understanding of the structure of the turbulence. The cross spectrum 
concerns itself primarily with the transport and conversion of energy 
and the transport of momentum and heat in the surface layer. All of 
these are fundamental in the understanding of micrometeorological 
processes in the atmosphere. The real part of the cross spectrum is 
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called cospectrum which for example, shows the effect of height above 
the ground on the scale of the eddies. The Imaginary part of the 
cross spectrum is called quadrature spectrum which gives some infor- 
mation about the vertical extent of eddies. Therefore the spectral 
behavior of atmospheric turbulence is of considerable and practical 
importance. 

Power and cross-spectral analysis of meteorological time series 
are generally based on the assumption that these series are stationary 
in the sense that their statistical properties are invariant with 
translations in time. At present, there exist techniques to analyse 
stationary time series records, but the techniques available for the 
analysis of non-stationary time series records are still inadequate 
and do not lend themselves to meaningful interpretations of physical 
problem. It is therefore necessary to use proper filtering techniques 
to adjust the non-stationary time series so that under certain cir- 
cumstances the existing statistical analysis for stationary time 
series may be used on these non-stationary data. It is also assumed 
that the time series are ergodic which permits time averaging to be 
used in place of ensemble averaging. 

Meteorological time series can mostly be represented as the sums 
of periodic (regular) and stationary (irregular) components since 
they almost always contain some definite periods, such as days, years 
or even decades etc. The longer than record length trends must be 
filtered out in order to eliminate the bias introduced In the spectral 
calculations. A moving-average and differencing high-pass filter Is 
applied to our time series so as to create a new, mean free, time series. 
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There are methods available to calculate the correlation functions 
via the fast Fourier transform method and to apply the thoroughly 
studied lag window to obtain the smoothed spectrum. For large volumes 
of data this is not considered to be an economical and practical 
approach. Development of direct Fourier transform methods as applied 
to meteorological data and consequent smoothing techniques are still 
in their infant stages and have not been thoroughly studied. The reason 
for this is partly due to the potential flexibility of the fast Fourier 
transform, which is much greater than the indirect or Blackman-Tukey 
method . 

The direct Fourier transform method or the periodogram approach 
is used for the calculations of power and cross spectra in this 
dissertation. A new algorithm is developed for the fast Fourier 
transformation which requires no bit-reversal procedures during the 
final stages of the calculations, A computer program based on non 
bit-reversal calculations is developed. Various time-domain smoothing 
functions have been applied to check the spectral variation. The 
calculated power and cross spectra agree fairly well with those re- 
ported in the literature of atmospheric turbulence. 

1.2 REVIEW OF LITERATURE 

The significance of correlation between the velocity of a fluid 
particle at one time and that of the same fluid particle at a later 
time or between simultaneous velocities at two fixed points was devel- 
oped in 1921 by G. I. Taylor [90], Taylor [91] further defined the 
scale of turbulence when it is applied either to the Lagrangian or to 
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the Eulerian concepts of fluid flow. The connection between the 
spectrum measured at a fixed point, and the correlation between 
simultaneous values of velocity measured at two points was also pro- 
posed by Taylor [92] In 1938. He pointed out that the correlation and 
the spectrum formed a Fourier transform pair. 

In general there exist two methods for the digital spectral 
analysis of stationary time series. The method which requires to 
calculate the correlation function first before taking its Fourier 
transformation to obtain the spectral density function is called the 
indirect method. Instead of calculating all the correlation functions, 
the direct method of computing spectral -density estimates is achieved 
by squaring the quantity which is obtained directly from the Fourier 
transformation of the raw observations. 

The spectral analysis using the indirect method has become a 
significant tool in the statistical analysis of stationary time series 
since the pioneer work of Daniel [22] and Bartlett [3,4]. They esti- 
mated consistent estimates of the spectral density function by mod- 
ification of the classical periodogram (i.e. harmonic) analysis. 

Optimum consistent estimates of the spectrum of a stationary time 
series were studied by Parzen [74]. Most of the work involving 
spectral analysis was done primarily by mathematical statisticians in 
the fifties. The practical situation of designing a spectral analysis 
satisfying many specific conditions had not been studied until Black- 
man and Tukey [10] in 1958. 



6 


With the advent of the electronic computers and improvement in 
numerical calculation techniques in spectral analysis, many valuable 
papers were published in the earlier Sixties (Jenkins[49], Parzen [75], 
Tukey [96] ). The users of the spectral analysis were still faced 
with the problem that if only a sample record of finite length is 
given, and no background information, is available, there still existed 
no precise method of obtaining an estimate which could be considered 
the most accurate one. One can actually construct many estimates of 
the spectrum by using different smoothing techniques, yet the optimum 
one, relies on the role of the chosen bandwidth [84]. 

The spectral calculations via correlation functions or Indirect 
method have been theoretically well established in the early Sixties. 

The periodogram approach or direct Fourier transform method in the 
time series analysis had been considered to be impractical because of 
the amount of computations involved. Actually, the periodogram as an 
estimate of the spectral density function of a stationary time series 
has had a long and controversial history starting with A. Schuster in 
1898 when he was investigating the hidden oscillations in meteorological 
phenomena with periods of 26 days. In 1965, Jones [53] reappraised 
the periodogram in spectral analysis and pointed out some advantages 
over the indirect method for multiple dimensional processes. 

During the controversy over the choice of direct or indirect 
approach to the calculation of the spectrum of stationary time series, 
the revolutionary paper on the Fourier transformation by Cooley and 
Tukey [19] was published in 1965. The algorithm they presented for 
the calculation of complex Fourier series permits the reduction of 
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operations by a factor of nV(N 1o 92 N) for a data sample length N. 

The latter should have a magnitude which is equal to an integer power 
of two. The operation consists of a complex multiplication followed 
by a complex addition. Numerous papers have been published to improve 
the algorithm or to explain the paper by Cooley and Tukey through 
different mathematical formulations. Ferrie [28] studied the relative 
advantages and disadvantages of various algorithms on the basis of 
execution time, storage and accuracy. He concluded that no single 
fast Fourier transform algorithm represents the best choice. 

Following the discovery of the FFT (Fast Fourier Transform) 
algorithm, earlier methods [10,75] no longer can be relied upon for 
the best statistical and computational procedures. Bingham et. al . 
[9], Cooley et. al, [18] and Welsh [105] presented various techniques 
for the estimation of power and cross- spectra via FFT, whose appli- 
cations in electronic engineering have been published in two exclusive 
issues in June 1967 and June 1969 of IEEE transactions on Audio and 
Electroacoustics. 

The direct and indirect methods do not produce the same results 
though both have relative advantages [27]. For a spectrum with a 
large slope, the direct method permits more window leakage than the 
indirect method. The indirect or correlation function method is more 
effective than the direct FFT or periodogram method in computing the 
spectrum of short time series. For large data samples, the direct or 
FFT method seems to be the only feasible one. The indirect method has 
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been considered not practical to handle a time series exceeding 
50,000 data points for spectral estimations. 

Oort and Taylor [69] applied the FFT technique to analyse the 
spectrum of the horizontal wind speed and to investigate the diurnal 
variability in the kinetic energy. Ryznar [85] measured wind and 
temperature profiles and turbulent fluctuations of wind velocity and 
temperature under various stability conditions of horizontally 
homogeneous turbulence in the atmospheric surface layer. He concluded 
that the integration .of the spectrum did not obtain the corresponding 
total variance prior to the spectral transformation. The reason is 
obviously due to the application of a "cosine bell" data window which 
reduced the variances of the spectral density estimate by a factor of 
3/8 [43]. Kaimal et. al . [55] calculated power spectra and cospectra 
of turbulence in the surface layer using the fast-Fourier transform 
technique. A new short time series was generated from the original 
time series by averaging the neighboring data points to investigate 
the information in the lower frequency range. 

Most of the books written on time series analysis such as 
Anderson [2], Grenander and Rosenblatt [41], Hannan [42], Wold [106] 
and the Brown symposium edited by Rosenblatt in 1963 [24] are con- 
cerned almost exclusively with theory. Practical applications of time 
series analysis started with Blackman and Tukey [10], though it was 
primarily from the point of view of the communication engineer. 

Books such as Brown [13], Fishman [31] and Granger and Hatanaka [39] 
provide a heuristic introduction to economic time series. The books 
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worthy of engineering applications are Bendat and Piersol [6], 

Jenkins and Watts [51], Otnes and Enochson [67] and the Wisconsin 
seminar edited by B. Harris [72] in 1966. 

The problem of obtaining the best spectral estimate via the fast 
Fourier transform lies in the selection of the smoothing function 
which depends also on the type of the data to be analyzed. In dealing 
with different types of data such as atmospheric wind data, biological 
data, acoustical data or economical data, the selection of the same 
smoothing function could lead to different interpretable results. 

Only by the proper choice of a suitable smoothing function in spectral 
analysis can meaningful results be achieved. 



CHAPTER II 


DATA PROCESSING SYSTEM 

The atmospheric wind turbulence data are collected through a 
three-dimensional split film anemometer (TSI-1080D) which consists 
of three split film sensors and a copper-constantan thermocouple 
for ambient temperature measurement. Each sensor consists of a 
0.006 inch diameter quartz rod coated with a platinum film of approxi- 
mately 1000 A in thickness. The total sensor length is 0.200 inches 
with approximately 0.08 inches of active usage. 

Data was sampled at a rate of 200 samples per second » and a 
special data acquisition and handling system has been developed to 
handle samples of half an hour duration. Detailed explanation of the 
data acquisition and handling system was presented in a report by 
Tieleman et, al . [95]. 

The data acquisition and data handling system was located in an 
instrumentation trailer positioned near the V.P.I. and S.U. low speed 
windtunnel. The instruments were used with a 350 ft. long connecting 
cable since mounting on the 300 ft. meteorological tower at NASA 
Wallops Station was proposed. After lengthy calibration in the 
low speed windtunnel, the probes were mounted in the atmosphere on top 
of the exchange section of this windtunnel. Since only a limited 
amount of time was available to test the probes in the atmosphere and 
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also since movement of the instrumentation trailer was a rather com- 
plicated operation, the instruments were tested in the atmosphere at 
the best location which could be reached with the 350 ft. cable from 
the position of the instrument trailer. The air flow on the top of 
the exchange section of the windtunnel (approximately 40 ft. above the 
ground) was a great deal influenced by surrounding buildings as well 
as the windtunnel and the exchange tower itself. 

In addition to the fast response T.S.I. probes, two slower re^ 
sponse Gill anemometers (see Appendix A) were positioned next to the 
T.S.I. probe so that comparisons could be made. The T.S.I. probes 
were positioned in a horizontal plane and could be rotated about a 
vertical axis by an antenna rotor into the prevailing wind direction. 

In the report by Tieleman et. al . [95] it was concluded among other 
things that the T.S.I. probes could operate with the best accuracy 

when they were directed into the mean wind. 

The two Gill anemometers were mounted also on the antenna rotor 
so that the three sensors could be moved all simultaneously. One of 
the Gill anemometers was mounted parallel to the T.S.I. probe and the 
second Gill anemometer was mounted perpendicular to the first one in 
the horizontal plane. Only two Gill anemometers were available so that 
only statistics from the two horizontal components could be compared. 


2.1 STATISTICAL ANALYSIS 



from the wind sensors were 


written on a 9-track 


magnetic tape in blocks with 209 samples. in each block. Computer 


programs were developed [95] to convert the raw data from analog 
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voltages to velocity components and temperature by using IBM 370/155 
digital computer. 

The computed velocity components and temperature were again stored 
on a 9-track magnetic tape with the block size increased twice to 418 
points in order that slight reduction in data-storage space on the tape 
and a saving of calculation time in the computer could be achieved. 

The data was collected at a rate of 200 samples per second, and 
therefore 360,000 samples were recorded during a half hour period. The 
sample size N is determined by the total recording time T and by the 
sampling Interval At through the following relationship 

N = ^ . (2.1.1) 

At 

The sampling rate, f., of the collected data is considered as an 
important factor in the data reduction process. This rate is the key 
to the determination of the Nyquist frequency, f„ which is the highest 
frequency one can obtain without introducing aliasing. The sampling 
rate and the Nyquist frequency are related by the following expression 



Let A.(n), 1 = 1,2,. ..,4 represent the three velocity components 
and the temperature recorded by the T.S.I. probe. The first three sub- 
scripts 1,2,3 denote the components of the instantaneous velocity vector 
in the coordinate system as determined by the directions of the sensors 
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on the T.S.I. probes (X, Y, Z). Subscript 4 represents the temperature 
at the point where the probe is located. 

The time series are divided into many blocks. The determination 
of the number of data points in each block depends primarily on the 
storage capacity of the available digital computer. The maximum number 
of data points should be chosen so that enough spectral information 
may be obtained in each block during subsequent spectral analysis. The 
basic requirement is that the selected block size should be long com- 
pared to the random fluctuations of the time history. - The sample mean 
and the sample probe yaw angle are calculated from the block mean and 
the block probe yaw angle. The mean value, the probe yaw angle, and 
the variance -in each block may be used to test the sample for 
stationarity as described in section 2.2. 


2.1.1 MEAN VALUES IN EACH BLOCK 

The sample is broken into M non-overlapping blocks each of which 
contains n points, so that Mn = N. The wind turbulence data may be 
considered to be a combination of a mean component and a fluctuating 
component. For a given number of n data points of the i-th time series 
Aj (n) in each block, the mean component may be described if ergodicity 
is assumed by the average value of all n values. 


1 n 


n 

j=l 


i=l....,4, 


(2.1.3) 


where subscript 1 refers to the number of the time series and superscript 
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m denotes the block number and the bar describes mean value. By 
choosing n = 8192 and M *= 44, we are analysizing 360.448 data points 
in approximately a sample of one half hour duration. The choice of 
8192 data points for each block is to meet the requirement for appli- 
ation of the fast Fourier transform algorithm that the number of data 
points to be transformed must be a power of two. Also storage re- 
strictions in the digital computer limits the block size. 

2.1,2 PROBE YAW ANGLE IN EACH BLOCK 

For each block of data points, a probe yaw angle is calcu- 
lated based on the previously obtained block mean values. The probe 
yaw angle is determined by the angle between the T.S.I, probes and the 
direction for which the lateral component of the mean velocity vanishes. 
This last direction fixes the so-called mean wind coordinate system. 

The x-direction the direction of the mean wind, the y-direction 
(lateral direction) perpendicular to the x-d1rection in the horizontal 
plane and the z-direction is vertically upward. 

The velocity components calculated initially in the coordinate 
system as determined by the directions of the sensors on the T.S.I. 
probes have to be transformed into components in the mean wind co- 
ordinate system for which the statistical quantities of the fluctu- 
ating velocities are of interest. 

The transformation of the velocity components from one coordinate 
system to the other is carried out in two steps. First, one obtains 
the components in the probe oriented coordinate system by the following 
relationship 
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S^*=E..a" (2.1.3) 

where 1=1 ,253 and j=l,2,3, are the direction cosines for the 
sensor oriented coordinate directions X Y Z and the probe oriented 
coordinate directions x*y*z*, the geometry of which is shown in Fig. 
(la). It is seen that the x* coordinate direction coincides with the 
probe axis, the origin of which is located at the probe tip. The y*- 
axis is perpendicular to the probe in a horizontal plane while z*- 
axis is vertically upward. 

Through geometric relations, Eq. (2.1.3) may be written in the 
following matrix form 


Al 


• 

0.57735 

0.57735 

0.57735 



* 

- 

0 

0.70711 

-0.70711 

1 


■jrfn* 

A 3 


-0.8165 

0,40824 

0.40824 



, 




. 


. . 


The second step is to transform to the mean wind coordinates. This 
further transformation between the probe oriented coordinate system 
x*y*z* and the coordinate system in the mean wind direction x, y, and 
z, as shown in Fig, (lb), may be achieved by the following relationship 

j=l,2,3 (2.1.5) 


or in matrix form 
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-m 

“i 


cosB 

-sinB 

0 


r * 1 

A"’ 

^1 

-m 

. 

sinB 

cosB 

0 


* 

^2 

-m 


0 

0 

1 


* 

»m 

^3 

k ^ Jl 


“ 






( 2 . 1 . 6 ) 


By combining the Eqs. (2.1.4) and (2.1.6) one can determine the 
X component of the block mean as follows 


= 0.57735 cosg, fl'j’ + {0.57735 cosg - 0.70711 sinB) 

+ (0.57735 COS6 + 0.70711 s1nB) S™ . (2.1.7) 

Maximizing the mean wind component in the x-d1rect1on with 
respect to the probe yaw angle (i.e., dO^/dB = 0, or similarly 
vanishing of U^) one obtains an expression for the block yaw angle 


tan 6 


-1,22457 



m=l ,2, . . . ,M 


( 2 . 1 . 8 ) 


The change of the probe yaw angle in each block indicates the 
shifting of the mean wind as averaged over 40.96 seconds (the 
latter is the total real time for each block). If there are signi- 
ficant changes in the block yaw angle from one block to the next, 
then we have reasons to believe the data to be nonstationary. The 
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same expressions 2.1.4, 2.1.6, and 2.1.8 may be used to determine 
the sample yaw angle and consequently the sample means, sample 
variances and sample covariances of the velocity components in the 
x,y,z coordinate directions. 

2.1.3 Block Variances of the Time Series 

The variances of the velocity components and temperature for 
each block are determined by the following equation 

^=1 " > A.(j) A (j) - . i=l,..,4 (2.1.9) 

where the subscripts 1,2,3 represent the three velocity components and 
subscript 4 denotes the temperature. The positive square root of the 
variance is called the standard deviation and is calculated by 

s'!' = {oi ) , 1=1, ..,4 (2.1.10) 

where m again denotes the block number. The block standard-deviations 
calculated by Eq. (2.1.10) will also be used to determine the non- 
stationarity of the data by using a statistical test which will be dis- 
cussed in section (2.2). 

2.1.4 Sample Mean and Sample Probe Yaw Angle 

The sample mean and the sample probe yavy angle are simply determined 
from the arithmetic averages of the block mean values and the block 
yaw angles in each sample by the following formulas 
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1 


S = - r A' 


m 


1=1, ...4 


( 2 . 1 . 11 ) 


and 

tanB = - 1.22475 (- ^ — ) , (2.1.12) 

where M is the total number of blocks or subseries. 

The sample mean value and the sample probe yaw-angle will be used 
to determine the velocity components in the mean wind coordinate 

t 

system by the transformation (2.1.4) and (2.1.6). 

Following the determination of the sample mean values for each of 
the three velocity components and the sample probe yaw angle in the 
data processing computer program called "DATPl", the trend and long 
period fluctuations will be removed through a program called "TREND". 

The theory of this filtering process will be described in section (2,3). 
The filtered data will be used to compute sample variances and co- 
variances as well as their respective values in the mean wind direc- 
tions through coordinate transformations. These computations will be 
included in another data processing computer program called "DATP2". 

2.2 Statistical Test 

The computed statistical properties for each block as described 
in the previous section 2,1 from a single sample record can be used 
to test the statlonarity of the time series. The wide variations of 
the probe yaw angle may indicate relatively large angular shifts In the 


wind direction . 
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If the statistical properties thus calculated do not vary signifi- 
cantly from one block to the next, one may be assured to a certain 
degree stationarity of the time series. Otherwise, the data of interest 
might contain nonstationarities and are not fit for further analysis. 

Without the detailed knowledge of the frequency composition of 
the calculated statistical properties, a non- parametric approach is 
necessary to determine whether or not the data are stationary [6]. A 
non-stationary trend test in either the mean or the variance is 
adapted from Kendall and Stuart [58]. The theoretical derivation is 
based on a paper written by Mann [65]. 

Suppose the block mean of each of the velocity components or the 
block yaw angle or their standard devtations are denoted by 

^1 * ^2 * ^3 ***** > 

where M denotes the total number of blocks in the sample. 

Now, a reverse arrangement of such a set of block variables Is 
defined as to occur every time 


Sj > , for all j > 1 

and i = 1 ,2,3, , M-1 . 

For a given value of index 1, the number of reverse arrangements 
for this given i is denoted by such that 

M 

Tt = E T- , (2.2.1) 

j=i+l 

where 

1 if Sj > S^- 


[0 if Sj < Si 
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Now, the total number of reverse arrangements Is expressed by 



^i 


M-1 M 
= E E 
1=1 j=i+l 



( 2 . 2 . 2 ) 


In the case the set of block variables shows an upward trend as 
far as their magnitude is concerned, the total number of reverse ar- 
rangements can be expected to be some relatively large number. Conversely, 
if the set shows a downward trend the total number of reverse arrangements 
is relatively small. If no discernible trend is present the total number 
of reverse arrangements is some intermediate number. 

The hypothesis Hq of no trend present against the alternative 
of a trend being present at an a-level of significance can be ex- 
pressed by 

Reject if T > H), T < t(l-^; M) 

(2.2.3) 

Accept if t(l- j; M) < T < t(^; M). 

In other words this hypothesis is similar to a symmetric 2-sided 
confidence Interval for T, with a confidence coefficient, 1 - a, 
and a lower limit, Ti_ = t(l-|; M). Therefore, for all T, 


Pj { Tl < T < T^j } 


1 - a . 


(2.2.4) 
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For a large number of block variables (i.e. larger than 10), a 
normal distribution for the total number of reverse arrangements of 
each set of block variables can be assumed. So that, the estimated 
value of T, T can be given by [64] 


T* . (2.2.5) 

where the average value and the variance of the total number of re- 
verse arrangements of each set are found to be 

. ( 2 . 2 . 6 ) 


and 


VAR[T] = 


2M^ + 3M^ - /5M 
72 


(2.2.7) 


respectively. These expressions are only valid when the set of block 
averages does not show a trend. 

The standard normal distribution function of the estimated values 
of T for an asymptotic N(0,1) distribution when is true, is given 

by 


^ [T*] = 


1 



( 2 . 2 . 8 ) 


— oo 
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where 


1 M(M-1 


* 

T = — r 


(2.2.9) 


+ 3H^ “ 5M 
72 


The right hand side of eq. (2.2.9) is seen to be in the standarized 
for. Of T* = Where s- = T . ^ 

If the normal variable S' is not "standard^.- its value must be sta,v 
darized according to T* = ^ ^ 5^^' 

The probability that the value of T is less than T for a given 
a-level of significance can be written as 


P[ T < T* ] = $ [T*] = a 


( 2 . 2 . 


where T* is the percentile of the standard Normal Distribution. By 
choosing different values of a, i.e. 0.95, 0.975, 0.99, the value 
can be obtained from the standard distribution table [TABLE I]. 

For different values of T*, a table of reverse arrangement distribution 
has been generated by using eq. (2.2.9). For different a-level of 
significance, and different number of blocks M, the values of 
t(^, M) and t(l - f; M) are given in Table II. 

If the value of the total number of the reverse arrangements falls 
outside the range in our criterion eq. (2.2.4) for a = 0.05, then a 
possible error of 5% can be made if the hypothesis that the data are 
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stationary is rejected. Clearly one wishes to commit such an error 
only rarely. If the value of T falls outside the range for a = 0.01, 
then one will have less chance to make this error, and consequently 
one has more confidence to reject the hypothesis that the data are 

stationary. 

This nonstationary trend test is generally effective in testing 
against linear or monotonic trends, as shown in Fig. 2 and ineffective 
against the type of nonstationarities as shown in Fig. 3 which show 
a reversal in the trend. The trend test is generally not successful 
In testing against the type of trends as shown in Fig. 4. 

2.3 Removal of the Trend 

Trend removal has been considered as an important step 1n the 
digital processing of random data. Large distortions can occur in 
the calculations of variances, covariances and spectral quantities if 

trends are not removed from the data. 

Trend removal is a special case of a general filtering process. 

Filters are designed to pass either low or high frequencies of the 
signal while attenuating or eliminating respectively high or low fre- 
quency components. Filters, which pass low frequency components are 
called low-pass filters while others which pass high frequencies are 
called high pass filters. A third type of filters which passes a 
band of intermediate frequencies and attenuate both very low and very 
high frequencies in the signal are called band pass filter. 

For meteorological time series, it is usually assumed that these 
series are statistically stationary and can be represented as sums of 
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regular (periodic) or irregular (stationary) components superimposed 
on each other. It is observed that these series almost always contain 
some definite periods that are due to certain external influences such 
as the time of day the data was taken and long periodic oscillations due 
to the presence of upstream obstacles. Since these periodic components 
may not enter into our statistical analysis of the observed time series 
some trend removal or filtering process is necessary. 

2.3.1 Methods of Trend Removal 

Different methods have been proposed for trend removal » the 
selection of an appropriate one depends largely on the practical sit- 
uations. For small samples, plotting is the best way to compare the 
filtering effects upon the unfiltered time series. For samples with 
a large number of data points, calculation time becomes one of the 
essential factors to make the selection. 

Dyer [25] suggested a modified difference filter applied to the 
meteorological data. However, the calculation of every auto-correlation 
coefficient between two neighboring values does not seem to be eco- 
nomical for the large volume of data encountered. Houbolt [47] pro- 
posed a high pass filter based on the idea of a symmetrical exponential 
filter. The requirement of making two passes in the calculation re- 
quires almost twice the computer time, which makes this method impractical. 

It is generally known that the best way of estimating a trend is 
to use a polynomial of low order by the least square method, yet this 
method does not represent the trend satisfactorily. As a result, one 
frequently makes use of the method of moving average based on the 
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idea that although a low-order polynomial, say a cubic polynomial, may 
not approximate the trend very satisfactory over the whole time Interval 
it may fit well over shorter intervals. 

2.3.2 High Pass Filter 

The high-pass filter applied to our data samples is obtained by 
first constructing the frequency response function of the equally- 
weighted moving-average low-pass filter which is the Fourier transform 
of the equally-weighted time function. The filter shape of the high- 
pass filter is obtained by squaring the frequency response function 
of the low-pass filter and then subtracted from unity. This particular 
high-pass filter is called the moving-average and differencing high- 
pass filter. 


2. 3. 2,1 VIeiqhtinq Function 

Digital filtering is simply a process by which a set of input data 
is transformed into a set of output data A° by means of a linear 
expression 





(2.3.1) 


where Wj^ are the suitably chosen weights and T is called the filter- 
ing interval. Eq. (2.3.1) can be regarded as the numerical approxima- 
tion of 


A°(t) = 


W(t) A^(t - t) dx 


(2.3.2) 


i , 

— CO 



26 


which states that the convolution of W(t) and A*(t) produces A°(t). 
The term dx is usually absorbed in the following form 


a° = WtA;t-^ ••• ^ W-1 "••• 

"2 




(2.3.3) 


where T must be an even integer. 

The weight W^, which is multiplied by the observation aJ, is 
termed the central principal weight. It is seen that the greatest 
weight Is placed on the most recent observation while both past and 
future observations receive syrnmetrically diminishing weights. In 
choosing an equally-weighted moving average filter* the weights selected 
are all equal to where T is the number of observations used in 
computing the mean or the filtering interval as described in Eq. (2.3.1). 
Therefore, the analytical form of the weighting function W(t) may be 
expressed by 


W(t) = 


r 1 

V 

0 . 



(2.3.4) 


It is seen that the weighting function W(t) Is applied to the 
observations from - ^ to ^ so it is an even function of t. The fre- 
quency response of the filter, H(f), is obtained from the Fourier 
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transformation of the weighting function W(t) by 


H(f) = 


W(t) e 


• 00 


(2.3.5) 


or because of symmetry 

oo 

H(f) = 2 f W(t) €05(2^^) dt. (2.3.6) 

Jo 

For a finite filtering interval T, Eq. (2.3.6) may be written as 
T/2 

H(f) = 2 I W(t) cos(27Tft) dt. (2.3.7) 

■^0 

Substituting Eq. (2.3.4) into the Eq. (2.3.7), the frequency 
response function of the moving average filter may be written as 


H(f) = 


simrft 

Tfft 


(2.3.8) 


which has a value of unity at f^^O. This frequency response function 
is calculated from the equally-weighted moving average time-function. 

To obtain a high pass filter, it is necessary to subtract the moving 
average value from the original data. 

The filtered data set for a moving average and differencing filter 
may be obtained by the following expression 
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T 

\ \ " "T+T \ ...,L-2 (2.3.9) 

i = -l 
2 

where is the filtered value at time t, is the original input 
value, T is the filtering interval and L is the total number of data 
points equal to twice the filtering interval. In this data processing 
system, I is the total number of the data in one block (i.e. 8192) 
while L = 16384. 

2. 3. 2. 2 Filter Shape for Moving-Average and Differencing High Pass 
Filter 

The filter shape for the moving-average and differencing high- 
pass filter is obtained from Eq. (2.3.8) by subtracting from unity 

1(f) = 1 - (2.3.10) 

which is shown in Fig. 5, 

The moving average and differencing filter obviously has the draw- 
back of failing to provide the trend values for the first half in the 
first block of the sample and for the last half in the last block of 
the sample. It Is not a great loss to have to forego the initial 
data values at the beginning of the first block, but the absence of 
trend values at the end of the last block is a serious handicap if we 
want to extrapolate into the future for forecasting. 



29 


2.4 Coordinate Traiisfonnation 

Four new, mean-removed time series a^, i=l...j4 have been created 
after having applied the moving-average and differencing high pass 
filter to the original raw observations. When i = 1,2 or 3, a. repre- 
sents the fluctuating velocity components in the sensor oriented co- 
ordinate system, and when i = 4, a^ denotes the fluctuating temperature. 

The sample probe yaw angle as calculated In Eq. (2.1.12) will be 
used in the calculation of coordinate transfonnation of the values of 
sample means, sample variances and covariances from the TSI or sensor 
oriented coordinate system to the mean wind direction xyz. The x- 
direction is the intersection of the vertical plane, which includes 
the total wind vector, with the horizontal plane. The y-direction is 
in the horizontal plane perpendicular to the x-direction and z-dlrec- 
tion is vertically upward, 

2.4.1 Variance and Covariance in each Block 

The filtered time series should have a near zero mean value, 
so variances and covariances in the blocked subseries m, m = 1 , 2 ,..., 

M-1 can be calculated by using the following definition 




= 1 
n 


n 

S a.(k) a.(k), i.j 
k=l ^ J 


1...4 


(2.4.1) 


where a^, i = 1,.,,4 represents the filtered time series. Equal subs- 
cripts i ® j in the above equation represent the variances while unequal 
subscripts 1 ^ j denote the covariances. The total number of blocks 
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have decreased by one due to the application of rnoving-average and diff- 
erencing high-pass filter. 

In order to correct the values of the block variances and co- 
variances of the time series, the filtered time series for each block 
should have exactly zero mean value. In this case, the following 
relationship should be used 


^ i 

1 j n 


n 

I 

k=l 


a^(k) a^(k) 




i,j=l,2,...,4 


(2.4.2) 


where the block mean value of the filtered time series is calculated by 



^ 2 a.(k) . 

^ k=l ’ 


i = 1,2, ...4 


(2.4.3) 


and m denotes the block number. 

Covariances of two time series measures the covariation between 
the related time series. The covariance will have zero value if the 
two time series are not related. 


2.4.2 Sample Mean Values 

The values of the sample means for the four time series have 
been calculated previously in Eq. (2.1.11) by using the unfiltered 
observations. The sample variances and covariances are calculated from 
the filtered series a., i = 1,2,,.., 4 by using arithmetic average 
of each blocked values as follows 

. 1 M-1 

‘irrr 


B 


1.j=l....4 


(2.4.4) 
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where superscripts m = 1,.., M-1 denotes the number of blocks to be 
averaged. 

The sample variances and covariances for the filtered time series 
a^. , i=l,..,4 with exactly zero mean value are calculated from 

a. a- * = J-i -- E aVa."'* . i,j=l,2,..4. (2.4.5) 

1 J M-1 1 J 


2.4.3 Mean Wind Direction Transformations 

The statistical quantities of interest should be expressed in the 
mean wind direction, therefore transformation of all the values of 
sample mean, sample variances and covariances are to be performed from 
the sensor oriented coordinate system to the mean wind coordinate 
system. Since the covariances between fluctuating temperature and 
three fluctuating velocity components are also important, a new trans- 
formation 4x4 matrix is obtained by combining the relations in Eqs. 
(2.1.4) and (2.1.6) as follows 


Ell 

^12 

El 3 

0 

^21 

hz 

^23 

0 

^31 

hz 

^33 

0 

0 

0 

0 

^44 


where = 0.57735 cos6 
= 0.57735 sin3 
E^3 = -0.8165 
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= 0.57735 cose - o. 70711 sin^ 

^22 ^ 0*57735 sine+ 0.70711 COS 3 
E 23 = 0.40824 

= 0.57735 cose+ 0.70711 sing 
E 32 = 0.57735 sine - 0.70711 cos$ 

E 33 = 0.40824 

^ 44 = 1 . 0 . 

The values of mean velocity components and temperature, variances 
and covariances In the mean wind direction can be obtained by the 
following transformations 

\ . i = (2.4.7) 

and 

“i*'j ° ®k®J, ■ i,j=l..--,4 (2.4.8) 

where the first three subscripts represent the three velocity com- 
ponents in the mean wind direction and the last subscript represents 
the temperature. The variance of temperature fluctuations in the 
sensor oriented coordinate system should be invariant under the 
coordinate transformations since temperature is a scalar quantity. 
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SPECTRAL ANALYSIS 

The estimates of spectral density functions and other spectral 
characteristics associated with stationary multiple time series are 
considered necessary in order to study the physical properties of the 
phenomenon in terms of its behavior in the frequency domain. The 
power spectrum shows how the variance or average power of the time 
series is distributed over the entire frequency range. The cross- 
spectrum describes the relationship between two time series in the 
frequency domain through determining the coherence function. Since 
the cross-spectral density function is a complex function, it can 
be expressed in terms of a real part (co-spectral density function) 
and an imaginary part {quadrature spectral density function) . 

3.1 Methods of Analysis 

In general, three different methods may be used to compute the 
power and the cross spectral densities. Each of the three methods is 
based on a different but asymptotically equivalent approach. These 
methods are 

a. the indirect or Blackman-Tukey method which takes the Fourier 
transform of the auto- or cross-correlation functions of the time 
series to obtain the required spectral density functions. 

b. the direct or fast Fourier transform method calculates a quantity 
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which is the direct Fourier transform of the time series. This 
quantity is then squared to obtain the spectral density functions, 
c. the filtering method can be used either on a digital or analog 
computer to obtain these spectral density functions. 

These three methods should turn out comparable results but these 
results are not necessarily identical, even if the same effective 
bandwidth at a given frequency is used. The three methods are suitable 
for computations of spectral density estimates, but they possess similar 
problems related to bandwidth, leakage, and statistical variability. 

3.1.1 Indirect or Blackman-Tukey Method 

This method requires the computation of the auto and the cross 
correlation functions before taking their Fourier transforms to obtain 
the spectral density functions. If the sampled observations {a^} and 
{bj} j = l,2,...,n come from two discrete time series with zero means, 
then auto- and cross-correlation functions can be calculated 
respectively as follows 


1 n-k 

= I' i ° (3.1.1) 

I 


and 


Cab(k) 


1 n-k 

n t t+k 


k > 0 . 


(3.1.2) 
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The power spectral density function of the discrete time series 
» j=l,2,...,n is estimated by 

J 

L-1 

G^(f) = 2At {C_(0) +2 I C ^{k)W(k)cos2TTfkAt}, 
a aa k~l 

0<f<^ (3.1.3) 

where W(k) is the lag window with truncation point L. The selection 
of an optimum lag window and best truncation point is usually done by 
trial and error procedures. 

The cross spectral density function for a stationary bivariate 
time series {a.} and {b.} , j=l,2,...,n consists of a real part of the 

J V 

cross spectral density function called cospectrum and an imaginary part 
of the cross-spectral density function called quadrature spectrum. 
Respective estimates of these spectra can be calculated as follows 


COgb(f) = 2At + 2 


L-1 

z 

k=l 


•tab(k)W(k) cos2TrfkAt} 




(3.1.4) 


and 


L-1 

= 4At E q3jj(k)w{k)sin 27rfkAt, 
k“l 




(3.1.5) 
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where i^ab(k) is the even part of the cross correlation function Cgb(k). 
That is 

Hab(k) = |■{Cab(k) + Cab(-k)), (3.1.6) 

and qab(k) is the odd part of the cross correlation function, namely 

“ 2 (3.1.7) 

The cospectrum describes the in-phase relationship of the two 
time series, while the quadrature spectrum depicts the out-of-phase 
relationship. The quadrature spectrum assumes zero value if the cross 
correlation function is even. The occurrence of a maximum correlation 
between the two times series {aj} and {bj}, at a non-zero lag will 
produce an odd function for Cg| 3 {k). 

The calculations of auto- and cross-correlations in equations 
3.1.1 and 3.1,2 involve a computational loop which is mainly a 
multiply-add operation, requiring execution time in modern high-speed 
digital computers of the order of seconds for just one single value 
for the time lag. In many physical applications, the number of time 
series to be analyzed can run into a fairly large number. Because of 
the time constraint, it is almost impossible to use this method to 
handle experimental records which run for half an hour and have a sample 
rate of 200 samples per second. As a result, alternative methods 
should be used for the spectral calculations of time series representing 
the ground level atmospheric turbulence. 
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3.1.2 Direct or Fast Fourier Transform Method 

Instead of calculating the correlation functions, this method uses 
the Fourier transformation of the discrete time series {aj}, j = 

1,2,... ,n directly. First one has to obtain the following quantity 


Hf = At i: f = 0,1,. ...n-1, (3.1.8) 

j=1 ^ 


where W(j) is the data window. Hence the power spectral estimate is 
obtained as follows 

^a(f) = Hf (3.1.9) 

* 

where is the complex conjugate of Hf. Equation (3.1.9) may also be 
written as 

®a(f) = sit • f = 0.1,2. ...|+ 1. (3.1.10) 

The introduction of the fast Fourier transform method makes the 
direct method extremely attractive in spectral density calculations. 

The detailed description of the application of the direct method in 
estimating power and cross spectral density functions is discussed in 
section 3.6. 


3.1.3 Bandwidth Filtering Method 

For a specific frequency Index k, it is necessary to assign 
frequencies {f|^} and bandwidth {Bj^l as follows 

° < fo < fl < ••• < fk ilk 

where At is the sampling Interval. The bandwidth {\} of a narrow band- 
pass filter may be described in three different ways, namely, the 
half-power point bandwidth, the noise bandwidth, or the equivalent 
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statistical bandwidth. For the half-power point bandwidth {B|^} used 
in this section, it is the frequency interval between the upper and 
lower frequencies where the filter attenuates an applied signal by 3 db 
below maximum transmissibility [6]. 

A separate bandpass filter, such as Chebyshev sine bandpass 


filter or Butterworth sine bandpass filter [67], is designed for each 

B. 

k. The filters have their half-power points at (f. - and 

Bk 

(f|^ + hertz. More specifically, the distance between the two half- 


power points is set to be hertz and f|^ is located midway between them. 

Let the output of the k^h filter be denoted by i=l,2,...,n 

which is filtered by using each of the filter k, then 




(3.1.11) 


where h^ and for j=l,2,...,t are the chosen filtering functions with 
total number of I weights. In practical numerical computations, it is 
noted that some input values and output values q . are necessarily 
set to zero initially. The power spectral density function is obtained 
by 



n 

z 

1=1 




(3.1.12) 


which means that the data are passed through a bandpass filter, squared, 
summed and finally normalized with the proper units. 
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3-2 Sam 0 l i hg Cons 1 derations for Random Data 

For the purpose of statistical analysis, most continuous sample 
records should be read at some fixed interval At and converted into 
digitized records for numerical calculations. Sampling defines the 
points at which the data are to be observed. Corresponding relation- 
ships exist between a random time sample record a(t) defined for the 
time interval from 0 to T seconds and its Fourier transform G(f) de- 
fined over a range of frequencies from 0 to F. However, both sample 
records a(t) and its Fourier transform G(f) are restricted by their 
respective time and frequency properties. Proper considerations should, 
therefore, be given to these problems in order to obtain better estimates 
of the spectral density functions. 

3.2.1 Resolution Difficulties 

Resolution is defined as the degree to which the true spectrum 
shows its narrow and tall peaks. Time and frequency are related inversely 
as can be seen from their physical units. The actual record lengths 
are finite instead of infinite in extent and the frequency bandwidth 
Af is also of finite width instead of near zero width. Due to this 
resolution problem, additional errors will be encountered in the estimates 
of the spectral density function. 

For a fixed record time T, the estimate of the power spectral density 
function might be improved by taking the frequency bandwidth Af by the 
following relationship 


T Af > C 


(3.2.1) 
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where C is a fixed constant. For a fixed record time T, values of 
Af which are too small will violate this uncertainty principle. 

One should attempt to increase T and at the same time decrease Af 
without corresponding increase in record length T will reduce the 
statistical reliability of the spectral density estimates. A decrease 
in bandwidth Af will improve the resolution but has to trade for 
statistical reliability. A compromise between a reasonable bandwidth 
and statistical reliability is, therefore, necessary. 


3.2.2 Aliasing Errors 


The question of the "aliasing" error arises as a result of 


sampling the data a(t) at equal Intervals of time At and later confusing 
some of the higher frequency contents in the original frequency space 
with the lower frequencies as can be seen in Fig. 6. The aliasing can 
easily be avoided electronically in the experimental system by filtering 
the signal before sampling so that the power above the maximum fre- 


quency fg is effectively removed. In digital data computations, care 
must be taken to avoid the occurrence of aliasing. 

If a data sample a(t) is sampled at equal intervals of time by 



(3.2.4) 


then fg, which is called Nyquist frequency, is the highest frequency at 
which spectral data can be attained without introducing aliasing 
errors. Any frequencies present in the data which are integer multiples 
of Nyquist frequency fg cannot be distinguished from fg. For the 
frequency f - nfg and At = it js seen that 
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sin2'FrfAt = sin2Tr(nfs)(^) = 0, for all n=l,2,... . (3.2.5) 

Similarly, the energy at an arbitrary frequency f cannot be 
separated from energies contributed by different frequencies such as 

f ± fg, f ± 2fg, f ± nfg, (3.2.6) 

and 

sin22TrfAt = sin22Tr(f + nfs)(-^), for all n=l ,2 (3.2.7) 

Thus, if frequencies higher than the Nyquist frequency fg are 
actually presented in the data, they will contribute their energies 
to lower frequencies with consequent errors in power spectral density 
estimates at these lower frequencies. 

To avoid this aliasing problem, one should choose the sampling 
frequency and sampling rate At in such a way that 

At = 2^1 2?^ — , (3.2.7) 

^'s ^‘max * 

where f^^g^ is the maximum frequency for which the data will be 
analyzed. It may also be written as 

fsiW- (3.2.8) 

It is, therefore, concluded that the time interval between 
successive samples should be such that the sampled data contain at 
least two samples per cycle of the highest frequency of interest. In 
case extraneous noise is present in the data samples, ten to twenty 
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samples per cycle of this highest frequency are considered practical 
unless the high frequency noise can be filtered off before the data 
are sampled. 

The effect of aliasing on the integral Fourier transform of a 
function a(t) can be shown as follows 


G(f) = 


. . -127rft^^ 

a(t)e dt 


(3.2.9) 


where i = /T and its inverse Fourier transform is 


a(t) = 




(3.2.10) 


The effect of sampling at finite intervals, evaluated at the 
points tj = jAt, j=0, ±1,±2, ... and F = ^ can be seen as 
follows 


a(tj) = 


i2-rrfj 

G(f)e ^ df 


Z 


G{f)e ^ df 
KF 


(3.2.11) 


The exponential function in the integrand of Eq. (3.2.11) is a 
periodic function of frequency f in the region from f = 0 to f = F, 
so it may be written as 


a(jAt) = 


■12TTjf 

Gp(f)e ^ df 


0 


(3.2.12) 
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where Gp(f) is the spectral density function with periodic function 
f such that 

CO 

Gn(f) = r 6(f + Kf) . (3.2.13) 

^ K=-^ 

The Gp(f) is different from the spectral density function G(f) 
defined in Eq. (3.2.9) because Gp(f) is the sum of the G(f)'s dis- 
placed by all multiples of F. This error is referred also as 
"aliasing" in the frequency domain. 

3.3 Tapering Function ■ Time Domain 

The selection of a tapering function is in many respects analogous 
to an engineering design of an electrical filter. Tapering is to 
multiply the time series by a data window, analogous to multiplying 
the correlation by a lag window in the indirect or Blackman-Tukey 
method. The method of using various lag windows to obtain smoothed 
spectral estimates have already been well established [77]. The 
problem of tapering the time series has not been thoroughly studied 
although different data window functions have been proposed [105]. 

By adopting the direct or FFT method for spectral estimation, 
emphasis will be placed on the use of data windows or tapering 
functions. The purpose of using a tapering function is to provide a 
slightly wider spectral window than would be obtained if a straight- 
forward harmonic analysis is used. From the time domain viewpoint, 
tapering is to round off potential discontinuities at each end of the 
time series. In the frequency domain, tapering is to suppress large 
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negative side lobes in the power spectrum. In general, tapering will 
change the mean and variances of the data sample since unequal weights 
are given to different portions of the time series. In effect, some 
data are lost and as a result degrees of freedom are lost. A scale 
factor is therefore necessary to compensate this difference in order 
to obtain the accurate spectral estimates.' 


3.3.1 Cosine Taper Data Window 

It is essential that the mean is removed from the data before the 
data window function is applied. Bingham et. al. [9] proposed a data 
smoothing function consisting of a short left-half cosine bell, a 
long constant and a short right-half cosine bell. Later they proposed 
a data window to taper both ends of the time series with the cosine 
bell each of which contains one tenth the time of the total sample 
time. The data between these cosine bells are multiplied by unity. 
This data window is called cosine taper data window and may be expres- 
sed in the following form 


Wt = (^) (1 - cos 


0 < t < O.IT 1 


Wt = 1 

Wt == (7) n “ cos 


O.IT < t < 0.9T 
0.9T < t < T 


(3.3.1) 


(or 0 < T-t < 0. IT 


where T is the total sample time. The corresponding smoothed filter 

shape Ct(f) which can be used in conjunction with FFT, is shown in 
2 

Fig. 8. This function is the Fourier transform of the data Window 

Wt, as is shown in Fig. 7. The Cj(f) has a wider main lobe with 

F 
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suppressed side lobes so as to reduce any possible power leakage and 
to prevent negative power estimates. 

The general form of the cosine-taper data window of Eq. (3.3.1) 
may be written as 


Wt = 1 (1 - cos 1^) 


(3.3.2) 

1 OTrt 

lOirt 


^t “ ^ ^ + 

e ^.)] 


lOnt 

20llt lOTTt 


Wt = e T - 

le T. 

(3.3.3) 


From the numerical coefficients inside the bracket of the above 
equation, it is seen that the cosine-bell data window is simply an 
extension of the Hanning weighting function with coefficients 
I* frequency domain smoothing. 

3.3.2 Other Data Windows 

There are few other proposed data windows [105] but they have not 
received wide attention. Their application in any practical problem 
has not been attempted. In applying these windows to the atmospheric 
turbulence data, the results predicted are generally lower in value 
as compared to those using cosine-taper data window which was discussed 
in the previous section. A brief discussion of these windows is 
presented in the following section. 
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3. 3. 2.1 Window No. 1 

Welsh [105] proposed a data window which is expressed by the fol- 
lowing function 


wi(t) = 1 - 


t “ 


n-1 




n+1 

2 


, t = 0,1,2,.. .,n-l. 


(3.3.5) 


The resulting spectral window corresponding to the data window 
(3.3.5) is given approximately by 


W](f) 


1 ^ 2 j- s1n(n+1 )Trf 

nU' ir2(n+l)f2 (n+l)Trf 


COs(n+l )'rrf]} 


(3.3.6) 


where 

1 n-1 2 

U' = ^ z w, (t) 
t=0 ‘ 


(3.3.7) 


and n is the total number of data points used in the computation and 
can be considered as a scale factor. To change n will mean a variation 
of the shape of the spectral window, W](f), by expansion or compression 
of the extent of the independent variable f. 

The half-power width is given by 


Aif = 


1.16 

n+1 


(3.3.8) 


When the same half-power width is used for comparison, the spectral 
window, W-|(f) turns to be almost identical to the spectral window as 
proposed by Hanning [10]. 
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3. 3. 3. 2 Mindow No, 2 

Welsh [105] also proposed another data window but it has a 
different shape as compared to window No. 1. 


W2(t) = 1 



t = 0,1,2, m,n-l . 


(3.3.9) 


The resultant spectral window W 2 (f) is again obtained by applying 
Fourier transform to the data window W 2 (t) and can be expressed as 
follows 


W2(f) - 



n+1 sin^{(n^-l )g^> 
^ {(n+1)|l}^ 


(3.3.10) 


where 

1 n-1 2 

U' = -1- i: w, (t) . (3.3.11) 

" t=0 

Changing the total length of the data set will again result in 
the change of the shape of the spectral window function W 2 (f) by 
expansion or compression of the extent of the independent variable f. 

The half-power width is found to be 

. - 1.28 (3.3.12) 

''2^ = IhT • 

When this half-power width is used to compare the shapes of the spectral 
windows, W 2 (f) is found to be very close to Parzen's spectral window 
[51], which possesses large negative side lobes. The presence of 
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negative side lobes in the spectral window function has great dis- 
advantages since prediction of negative power spectral density function 
may occur. 

3.4 Smoothing Function - Frequency Domain 

The smoothing of the spectral estimates in the frequency domain 
may be achieved by either averaging the estimates at the corresponding 
frequency points of segmented subseries in a time series (segment 
averaging) or averaging the spectral estimates among neighboring points 
in the spectrum function (frequency smoothing). The combination of 
both segment averaging and frequency smoothing may be applied to obtain 
a smooth spectral estimate. This type of smoothing is referred to as 
combined averaging. 

A plot of the individual power estimates versus frequency for a 
one half hour sample will be next to impossible. Even the plotting 
of spectrum functions which are smoothed using the segment averaging 
method will show many small individual peaks. These small peaks and 
valleys are insignificant in the explanation of the turbulence 
structure, because they may represent sampling fluctuations rather 
than any systematic physical variations. Frequency smoothing may 
average out these peaks in order to obtain a more useful representation 
of the spectrum but has limitations as far as the resolution is con- 
cerned. However, for a long time series covering nearly four decades 
on the frequency scale, the resolution is generally two or more orders 
of magnitude greater than actually required [69], 
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3.4.1 Segment Averaging 

The data sample of one half hour duration has been previously 
blocked into M different subseries each of which has n data points. 

The total number of blocks has been reduced by one as a result of 
applying the moving-average and differencing high-pass filter. By 
using the fast Fourier transformation technique, there will be M-1 
number of blocked spectral estimates each of them of length n. The 
M-1 spectral estimates are averaged over corresponding frequencies 
to obtain a smoothed spectral estimate given by 

, M-1 

G(f^) = %(fi) 1-1,2 n. (3.4.1) 

The application of segment-average smoothing will increase the 
effective resolution bandwidth Be depending on the number of blocks 
to be averaged. The spectral window before applying the segment 
averaging was triangular in shape with Be = ^ . After applying the 
segment averaging, the spectral window Is stHl triangular in shape 
except the effective resolution bandwidth Be will be wider (i.e. Be = 

Ml 

as shown in Fig, 9. The spectral estimate 6(fi) In Eq. (3.4.1) 
may be considered as representing the midpoint of the frequency inter- 
val covered by Bg. A total of n spectral estimates can be obtained. 

3.4.2 Frequency Smoothing 

After averaging the i neighboring spectral estimates of the power 
spectrum, -the new spectral function can be expressed in terms of the 
original one as follows 
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G(f,<) = i s G(fi) K = 1,2,... " 

^ ^ 

where £ is the number of neighboring frequency components to be aver- 
aged and n is the total number of original spectral estimates. The 
spectral window before frequency smoothing was triangular in shape with 
the effective resolution bandwidth Bg = y , After frequency smoothing, 
the spectral window will be trapezoidal in shape with much wider ef- 
fective resolution bandwidth (I.e. Be = j) as shown in Fig. 10. The 
spectral estimate may be considered as representing the midpoint of the 
frequency interval from f(( to f«+jt*i» where K denotes interval number. 

3.4.3 Combined Averaging 

The smoothing of spectral estimates can be achieved more effective- 
ly by first applying segment averaging followed by frequency smoothing 
which is known as combined averaging. As a result of applying the 
combined averaging technique for smoothing of the spectral density esti- 
mates, the final effective resolution bandwidth becomes much wider and 
can be approximated by 

Be =T^ (3.4.3) 

where M 1s the number oSf blocks to be averaged and z is the number of 
averaged spectral estimates and T is the total sample time. The number 
of degrees of freedom y is 

Y = 2BeT = 2M£. (3.4.4) 

which can be interpreted as the total number of real and imaginary 
components within the bandwidth, Bg. Since each of the spectral esti- 
mates is in itself a Gaussian random variable, the squaring and 
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adding of them will produce a Chi-square distribution. The broader 
effective resolution bandwidth gives the reduced normalized standard 
error which will be explained in the next chapter. 

3.4.3 Proposed Frequency Smoothing Technique 

The frequency smoothing is usually done by choosing an equal 
number of neighboring frequency components of the spectral estimates 
to be averaged. This technique is effective only if the power spectral 
density estimates are distributed evenly throughout the whole frequency 
range of interest. If the power spectral density estimates are con- 
centrated either in the low or in the high frequency range, a new 
frequency smoothing technique is proposed in order that more informa- 
tion is to be obtained In the range of Interest. 

In order to obtain more spectral information in the lowest 
frequency range of the atmospheric turbulence spectrum, it is proposed 
that the total number of spectral estimates after segment averaging 
be separated into different averaging regions. 

Since the power spectral estimate at zero frequency is of no 
significance for reliable analysis, the smoothing starts with the 
second value of the spectral estimates after zero frequency. 

The power spectral density estimates at the lowest frequency 
range have been obtained by only going through segment averaging 
(i .e. = 1 , n = 8 in Eq, (3.4,2)). The average value i used in the 

frequency smoothing method is chosen in a manner of geometric progres- 
sion (i.e. i = 2^,2^,2^,,.,.). For different values of i different 
spectral estimates will be obtained as expressed by the following 
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G(fk) =1 j G^(fk) 

X . 1 = 1 

= 4. 

G(fk) = I X Gi(fk) 
1 = 1 

Jl = 16 , 

G(fk) =I Gi{f^) 

1=1 

Z = 64, 

G(fk) -J 6l(fk) 

and z = 256, 

G(fk) = T s,-(<'k) 


k = 1,2,... ,8 . (3.4.5) 

k = 9,10,... ,38 , (3.4.6) 

k = 39,40 62 , (3.4.7) 

\ 

k = 63,64,..., 86 , (3.4.8) 

k = 87,88 94 . (3.4.9) 


In this case, 4097 unique values of spectral density estimates Gi*(fj() 
have been used to obtain 94 final smoothed spectral density estimates 

A 

6(f(^). The overall filter spacing for this proposed frequency smoothing 
technique Is shown in Figure 11. 

3.5 Fast Fourier Transform 

The fast Fourier transform (FFT) is an efficient and time saving 
algorithm for the computation of the harmonic amplitudes in the 
frequency domain from uniformly spaced input data points. It is used 
to analyse the periodic phenomena of a time series as it converts to 
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a frequency function. The time savings of this algorithm will permit 
the reduction of operations by a factor of N^/N 1093 ^ number 

of data points in the sample is a power of two. 

Many widely used FFT algorithms involved bit-reversal procedures 
to sort the Fourier coefficients into the proper order in the final 
stages of the transformation. The bit-reversal procedures can also 
be avoided if the input data have been sorted into different orders 
so the output coefficients will be in the proper order. 

Uhrich [99] developed a Fortran program without requiring bit- 
reversal procedures but lacked a supporting theory. His algorithm 
requires a 2xN array for variable storage in the computer rather than 
a IxN storage for bit-reversal sorting. The required computation time, 
therefore, is almost twice as long as Is necessary. 

A different algorithm for FFT calculations without bit-reversal 
and sorting procedures together with complete mathematical formu- 
lations will be presented in the following sections. 

3.5.1 Properties of Fourier Transform 

In order to understand the fantastic time savings involved in 
FFT computations, some basic properties of the Fourier transformation 
need to be explained. In processing digitized signals or samples, 
the attention will be focused on the properties of the discrete finite 
Fourier transform (DFFT). Some useful properties related to the 
development of the fast Fourier transform will also be discussed. 

Property 1 : The Fourier transform transfers the N data points of 
the time series (aQsa-| ja^,. to the N-spectral values 
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,H2». . . ) in the frequency domain, as defined by 

N-1 

Hy.(f) = E a^(t) r = 0,1,2,... ,N-1 (3.5.1) 

n=0 

-i27r 

where W = e ^ • A very similar equation defining the inverse finite 
Fourier transform, which transforms the discrete spectral values back 
to the original discrete time series can be written as 

1 

an(t) = ^ S H-(f) n = 0,1,2 N-1. (3.5.2) 

r=0 

In case the digitally recorded data points a.,a. . ,a,, , are 

0 1 N-1 

finite but non-periodic (random), the sample values as defined by Eq. 
(3.5.2) can be proven to be periodic and infinite as follows 

r=0 

or 

^ 2 H^(f) W (W*^) (3.5.3) 

r=0 

where W is the root of 1 or = 1 , so (W^) = 1 and it follows 

that 

N-1 

an+N(t) = f E H^(f) W'*"” = a^(t) . (3.5.4) 

r=0 

Similarly one can prove that 
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Property 2 : If the time series (aQ.a-j ) is real, then 

the second part of its finite Fourier transform is equal to the 
complex conjugate of the first part. The spectral values Hq.H-j,..,, 
can be separated into two parts 


, Hfj, Hn .Hu Hf,_i 

2 r"’ 2^2 


First Part Second Part 

It is necessary to prove that 


Hj^_j(f) = Hj(f) for 0 < j < N-1, (3.5.5) 


where the star superscript refers to the complex conjugate. From Eq. 
(3.5.1 ) , one can write 


Hf^.j(f) 

Using the property 1 , 
HN-j(f) = 


'z’ a„(t) W^N-j)n . 
n=0 

Eq. (3.5.6) reduces to 

r a„(t) W’J" . 
n==0 


Now, 




1 


COSI 


2ti 


. 2tt 

slnjj- 


(3.5.6) 


(3.5.7) 


(3.5.8) 


or 

,,-l 2TT . . . 2 it 

W = cos]^ + 1 sinjp = W . 


(3.5.9) 
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Hence Eq. (3.5.6) can be written as 

* .*n * 

= s a (t)(W^”) . (3.5.10) 

n=0 

For a real valued time series* a^Ct) = a*(t) and Eq. (3.5.10) can be 
written as 


HN-j(f) = 


N-1 

s: a„{t) WJ" 
n=0 


H.(f) . 

J 


(3.5.11) 


Property 3 : The atmospheric turbulence data are always a sequence 
of real numbers. By applying Fourier transformation we usually store 
the actual (real) input in the real -part array of the computer and 
zeros in the imaginary-part array of the computer. This requires 
storage length of 2H and produces 2N components of N complex Fourier 
coefficients. In fact, we used only N Fourier coefficients which 
should take much less storage in the computer. By forming a complex 
time series (Cq,c-| , . . . ,Cnj_i ) from two real input time series 
(aQ,a^ ,. . . ) and (Bq, 6^ , . . . , one can perform Fourier trans- 

formation utilizing complex input in a form of = a^, + i for 
r - 0,1 ,2, . . . ,N-1 and i = . Then the combined spectral series 

will be in the form of 


Cr = \ + i By. r = 0,1 .2, . . . ,N-1 (3.5.12) 

where Ay, and By. are both complex. Let 
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and 


= a-, + 1 a« 
‘r 


By» — b*j + i bo 

'r > 


(3.5.13) 


where a-i ,ap »b-i , and b« are real numbers, 
r r r ^r 

By using property 2, one can write 


^N-r = 1 


{3.5.14) 


where and By. are complex conjugate of Ay, and By, respectively. 
Both Eqs. (3.5.12) and (3.5.14) can be written in form of 


and 


Cr “ (fl] ^ 32 ) ^ (^1 ^ b ? ) 

r r r r 


^N-r " ^ '* ■ '' ^2^) • 


Hence it follows that 


and 


Cy, - (ai - bp ) + i (ap + b-i ) 

r r r 


^N-r "" + bp ) + i (ap - L ) 

'r ‘r 


(3.5.15) 


(3.5.16) 


By combining Eqs. (3.5.15) and (3.5.16) one obtains 


Cy. + Cw_y, ~ ~ 


(3.5.17) 
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Cr - 


21 




b2^ = Br 


(3.5.18) 


where and B^* for r = 0,1 ,2,. . - ,N-1 are the Fourier transform 
values of the two real time series a^ and bp respectively. Both Eqs. 
(3.5.17) and (3.5.18) can be written in more easily perceptive forms 
as follows 


Ap(f) = Re 


(Cr(f) + CN_^(f))' 


(Mf) - 

L 2 J 

+ 1 iin 

CO 
1 


(3.5.19) 


and 



(Cp(f) + CN.r(f)) 


(Cr(f) - CN.r(f)) 

CM 

- i Re 

1 

CM 


(3.5.20) 


There are many other properties (i.e. circular convolution of 
two time series) which may be found in most of the standard textbooks 
on Fourier transformation. 

3.5.2 Basic Theory of Fast Fourier Transform (FFT) 

The basic theory of FFT was developed by Cooley and Tukey [19] in 
a subtle way so that it was somewhat difficult to understand. The 
main idea involved in the entire development was to continuously 
reduce a time series into two final point functions. It may also be 
illustrated by the principle of matrix factorization [34]. The theory 
presented here is to illustrate in a simple mathematical form the 
reasons for the time-saving capacities of this algorithm. 
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In the sampled-data case, the discrete Fourier transform (DFT) 
is defined previously by 

N-1 

= i: r = 0,1 ,2,. . . ,N-1 (3.5.21) 

k=0 

where Hy- is the coefficient of the DFT and a|( denotes the 
sample of the time series which consists of N samples and i = /T. 

The definition of DFT is not consistent throughout the literature. 

Some authors prefer to use H^"/N as DFT coefficients, others use 
Therefore, care should be exercised in doing the numerical 
computations. The ai^'s in expression 3.5.21 can be either real or 
complex but the H^'s are almost always complex. Eq. (3.5.21) may be 
written in a simplified form 

rk 

Hy, = r r = 0,l,2 N-1 (3.5.22) 

k=0 

where or W = , 

3.5.2. 1 Illustrated Example 

To become more familiar with the discrete fast Fourier transfor- 
mation as defined in Eq. (3.5.22), an example will be used to describe 
this technique for 8 data points. For a random digitized time series 
a|c* k = 0,1 ,2,..., 7 the Fourier transform of Eq. (3.5.22) may be 
written as follows 
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^0 ~ ^1 ^ ^2 ^3 ^4 ^5 ^ ^ ^7 

H-] = ag + a-jW + + 331*1^ + a^W^ + a5W® + agW® + afA^ 

H2 = ag + aiW2 + a2W^ + a3W6 + a4W + a5W2 + agW'* + ayW® 

H3 = ag + a^H^ + a2W6 + a3W^ + 341/ + asW^ + agW^ + a7W^ 

f 3 5 23 ) 

H4 = aQ + a^W^ + 32 + S3W^ + 34 + a5W^ ■*■ 3g + ayW^ 

H5 = ao + aiW^ + aaW^ + 33^!^ + a4W^ + agW^ + aeW® + a7W^ 

^6 ” ®0 '*' ^1^^ 32W^ + a3W^ + 34 + agW^ + agW^ + a7W^ 

H7 = ag + a-]W^ + 32W® + a3W^ + a4W^ + agW^ + agW^ + a7W^ 

where use has been made of the relations = 1 , = W*^ and = -1. 

It Is clear that the complete calculations of Eq. (3.5.23) requires 64 
complex multiplications and additions. Samples which consist of 
thousands of data points would require an extremely large number of 
computer computations and also require a large storage. But from Eq. 
(3.5.23) some basic characteristics of the Fourier transform may be 

1 

recognized. Obviously, symmetry in the right hand side of theie 
expressions can be observed. The first equation in (3.5.23) is simply 
the sum of all the data points. 

The FFT method Is essentially to divide the total number of data 

points in half which gives two sequences and then dividing these 

sequences in half again to give four short sequences each consisting 

of two terms. It will be shown that the shorter sequences require 

fewer operations than the longer sequences. , It Is profitable to 

separate the original data points into two shorter sequences bj^, 

N 

X, = 0,1 ,2,. composed of only the even-numbered data points, and 
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s, = 0,1,2,..,,^ composed of only the odd-numbered samples. Now 
the discrete Fourier transform of the shorter sequences for the 8-data 
points sample is given by 

3 2rk 

By. = z r = 0,l,2,3 

k=0 

and (3.5.24) 

^ 2rk 

Cy, = I CkW r = 0,1.2,3, 

k=0 

where = cos(|^) - i sln(^). The computing time for obtaining the 
Fourier coefficients and Cy,, r = 0,1, 2, 3 Is now reduced to 2 ( 4 )^ = 

32. By expanding Eq. (3.5.24) one should see more clearly the ad- 
vantages of separation of the long sequence, 


Bq = bg + b^ + b2 + = 9 q + 83 + + ag 

Bi = bQ + b-]W^ + b2W^ + b^W^ 

^2 “ ^0 b^W^ 

B3 = bo + b-iW^ + b2W^ + b^W^ 

Cq = Cq + Cl + C2 + C3 = + 33 + ag + aj 

Cl = Cq + c^w^ + C21/ + C3M® 

C 2 = Cq + ciw^ + C 2 + C 3 / 

C3 = Cq + C-jW® + CgW^ + c^w^ . 


(3.5.25) 


To illustrate the advantages of FFT more clearly, the spectral 
values k = 0,1 ,2,. . . ,N-1 can be expressed in terms of spectral 

M 

values B(|^ and Cj^, i - 0,1,...,^ as follows 
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Hq = Bq + Cq 
H] = B-] + WC-| 

H2 = B2 + V|2C2 
H3 = 83 + w3c3 
H4 = Bq + W^Cq 
H g = Bi + 

Hg = B 2 + W®C2 

Hy = 83 + W^C3 . 


(3.5.26) 


Now the sequences b£ and C£, i - 0,1 ,2,3 can be further halved to 
obtain four short sequences of two terms each. They may be written as 
fol 1 ows 


1 

Dy. = I d4W 
1 

Ey, = 2 e^W 

0=0 ^ 

1 , 

Fr = I fiW 
0=0 ^ 


,4rj 


4rj 




r * 0,1 


r = 0,1 


r * 0,1 


(3.5,27) 


I 1 ' 4ri 

^ j=0 ^ 


r = 0,1 


By expanding Eq. (3.5.27), they are given as follows 
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Dq = do + di = Sq 

D-j = do + d-|W^ = 3 q + 

Eq = eo + ej = a2 + ag 

= Sn + eiW^ = ao + 

, \ ^ ^ (3.5.28) 

Fq = fo + f] = ai + ag 

Fi = fo + = ai + agW"* 

Kq = ko + = 83 + a; 

K-j = kQ + k-|W^ = 83 + a-^W^ 

Combining Eqs. (3.5.25) and (3.5.28), the spectral values By. and 
Cy,, r = 0,1, 2, 3 may be expressed In terms of the four short spectral 
sequences as follows 

Bq = ^0 

Bi = 

82 = Do + EqW'^ 

Br, = Di + EiW® 

, (3.5.29) 

Cq = Fq + 1^0 

Cl = F-i + Ki'w^ 

C2 = Fg + KqW^ 

C3 = Fi + Ki'w6 

where relations W® = 1 and are used. Therefore, the basis to 

calculate the finite Fourier transform of a time series aj^, k = 0,1,2, 
...,7 is to use the sets of equations (3,5.28), (3.5.29) and (3.5.26) 
in this order. 
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The example Illustrated above gives the basic reduction procedures 
needed to calculate the finite Fourier transform of a time series with 
eight data samples. From Eqs. (3.5.26), (3.5.28) and (3.5.29), each 
one of them requires 8 complex multiplications and additions. Thus 
the finite Fourier transform calculations can be achieved in a total of 
24 complex multiplications and additions. Using the straight forward 
calculations in Eq. (3.5.23), a total of 64 complex multiplications 
and additions are required. 

The reduction of the time sequence can be carried on as long as 
each reduced data sequence has a number of data points that is divisible 
by two. This FFT technique can be generalized to handle a time series 
of any length N as long as N = 2^, where p is any integer larger 
than 1. The number of complex multiplications and additions is pN 
instead of and therefore since always p <5 N the number of calcu- 
lations is drastically reduced. 

In suirenary, the basic steps for using FFT algorithm to perform 
spectral calculations for a time series with N = 8 may be written as 
follows; 

1. Use Eq. (3.5.28) to calculate the spectral values (Dq, D] , Eq, 

I I 

E] > ^1> Kq, K]). These eight values may be stored in the position 

(complex) occupied by the original data values for they are no longer 
needed. 

2. Use Eq. (3.5.29) to calculate the new set of spectral values 

(Bqs B-j , 82® B3, Cq, C-| , C2, C3). 
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3. Use Eq. (3.5.26) to obtain the required finite Fourier 
transform set of spectral values (Hq, H4, K2, Hg, , H 5 . H 3 , Hy). 
It is seen that these spectral values are not arranged in the 
required order. A bit-reversal procedure is therefore necessary to 
align them in the proper order. 


3.5. 2. 2 Mathematical Formulations 

There is a large number of fast Fourier transform algorithms and 
computer programs available, the details of their computations are 
difficult to understand. It is therefore advantageous to develop a 
new algorithm based on the original theory presented by Cooley and 
Tukey [19]. The general mathematical formulations for FFT without bit- 
reversal procedure are presented in this section. A computer program 
based upon the mathematical formulations of FFT is developed and pre- 
sented in the appendix. 

From the definition of the discrete Fourier transform in Eq. 
(3.5.22), r may be expressed in the form of 


r = 


^p-i^ 




J'o 


(3.5.30) 


where = 0 or 1 for all t. 

The time series a|^, k = 0,1 ,2, . . . ,N-1 can be separated into two 
time sequences each containing half of the total number of data points, 
thus Eq. (3.5.22) becomes 

Ni-1 ^ N-1 ^ 

= E aj.Wn^ + S auWn” 
k=0 ^ k=Ni " ® 


(3.5.31) 
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where Wq = e i = N-| ="2 = 2 ^^ and p is the total number of 

times necessary to reduce the time series into two-point functions. 

By changing the sunmation limits in the second term on the right hand 
side of Eq. (3.5.31), it may be written as 

N 1 -I r 

= s [a,^ + . (3.5.32) 

0 

Furthermore^ Eq. (3.5.30) can be rearranged as follows 

— 2"^ “ ■^p-1^^ ^ + ... + 3*22 + j-| E r-j . (3.5.33) 


Substituting the newly defined value r^ in Eq. (3.5.33) into the Eq. 
(3.5.32) one obtains 


H 


r 



k=0 


r ^ 

M 

2 

■ ®k + ak+Ni'^0 

m \ 

• Wq ^ 


r-|k 

Wo^ » 


(3.5.34) 


where use has been made of the fact that = 1 for the integers 
n = l,2,...,p. Eq. (3.5.32) may again be written as 


r-jk 

H Wn^ k=0,l,2 N,-l (3.5.35) 

" k=0 


where 


(1) 

®2k+Jo 




(3.5.36) 
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By using similar procedures* Eq. (3.5.35) can be further reduced 
into another set of time sequences as follows 

r-| k rn k 

, = r a-. . U- 

r 


or 


H. 


N 2-1 0) '^5 ^7 

k=0 “ |c=n2 0 

N 2 -I (1) Jl^ N2-1 (1) 

^ k^O “2k4o“o ’ " 

J'll 


Nl 


Nj-I 
Hr= Z 
k=0 


(1) ^ (1) „ 2 

®2k+jg ®2k+jo+Ni“o 


^ Wo 


j-j k 

NT 


r k 
r^K 


Wo 


where N 


_ W] . P-2 N 
2-^-2 - 4 


'•2 


*^1-Jl . .sP-3 . p-4 

= ^ ^p-1^ ■** ^p-2^ + ... + + j ‘2 


j 


'*0 = (-1)^^ . 


{3.5.37) 


Eq. (3.5.37) can again be written as 


T2k 

■ k^o “4k.20i+jo“'>”^ 


H - 

Hr " 2 S 


k = 0,1,2,... ,N2~1 (3.5.38) 


where 


( 2 ) 


''''' = 

4k+2ji+jo ^ 2k+j(j * ®2k+Jo+N/‘'^ ^Wg 


j^k 

NT 


(3.5.39) 


By employing similar halfing procedures, a general form may be 
developed through induction as follows 



y 
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Hy, - E abWrj 

k=0 ^ 

IiJi 

^1-1 (1) , Ni 

= k=o "''"Jo 0 

r^k 

_ (2) Nf 

= ^Eq ®4k+2ji+j()W0 

r^k 

■ k^o ®8k+4j2+2ji+j(,Wo 

Y' k 

, Vri (p-i) ffr 

k=0 ^2P'’k+2P'2jp.2+2P‘3jp.3+...+2jT+jo”° 

(3.5.40) 

where p H the number of times the data sequences are decomposed into 
shorter sequences. 

By employing halting procedures a total of p times, one obtains 
the following expression 

. V’ 

^ k=0 *2Pk+2P-ljp.i+2P'2jp.2+...+2j,+jQ'^0 P , 

k = 0,1,2,. ...Np-l (3.5.41) 


where 
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(P) 

2Pk+2P-''j ,+2P-2j 2+. 


. .+2j]+jo 


1 

r 

q=0 




2P‘^k+2P-2jp.2+2P-2jp.3+...2ji+jo+qNp.i 


■^p-1 


N 


P 



The summation limit in Eq. (3.5.41) is no longer in existence if the 
time series of 2^ data samples is halved continuously a total of p 
times. Thus, Eq. (3.5.41) can be written as 


H 


r 


(p) 

2^ ^Jp-T+2P"^jp.2+..*+2ji+jo 



(3.5.42) 


which is the final result of a single point transformation. 

The previous mathematical developments may also be written in 
another form of notation so as to be able to grasp more clearly the 
ideas involved in developing the FFT computer program. In general, 
Eq. (3.5.22) may be written in the following form 


N-1 

H(r) = i a(k , 
k»0 P"' 



rk 

• . . tk-j >kQ)WQ 


r = 0,1,2,...,N-1 


(3.5.43) 


where r and k is expressed by 
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+ ••• + 2Ji + jo • 

(3.5.44) 

k = 2P-lkp.i + 2P-2kp_2 + ... + 2ki + ko 

and kj^ and jj^ take on values for 0 and 1 only, for all fi.. 

Expanding the exponential term on the right hand side of Eq. (3,5.43), 
it gives 


^0 


rk 

N 




r2P-''kp.i/N y rk^^VN 


rk(l) 

= Wg'^oVl''? WgN 

= Wg'’lVl WgnK^’VNi (3.5.45) 


where 


and 


N, =^2^ . 


r-jg 


ri = — = 2P-2jp.., + 2P-3jp.2 + ... + 202 + Ji 


k(l) = 2 P- 2 kp _2 + 2P-3kp_3 + ... + 2ki + kg 


(3.5.46) 


Equation (3.5.43) may 

1 1 

H(r) = I z 

kg-O k-|=0 


be expanded into the following form 
1 1 


I 

k2=0 



a{kp_i ,kp_2>* 


* *^1 


Wgjo''p-l''2 „^nk(l)/Ni (3.5.47) 

where Wg'“'''‘P-1 1s unity due to the fact that both ri and kp_, assume 
an Integer value. 
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Equation (3.5.47) may be written as follows 

®1 ^*^p-2**^p-3 


1 1 1 

H(r) = 1 l l 

kQ=0 k-|=0 k2=0 

rik(l ) 




Wo 


N 


(3.5.48) 


where (kp_ 2 .kp, 3 .. . ■ ,k^ ,ko;3*o) 


1 . jol<D-l/2 

2: a(kn-7 ^ '^O 

Vr° 


(3.5.49) 


Similarly, by letting 


^2f*^p-3*^p-4’* • * **^1 *^0^ 

= ^ ai(kp-2,kp.3,...,k^,ko,jo)Wo Wq (3.5.50) 

Eq. (3.5.47) may be further reduced to 


1 1 

H(r) = E 2 

kQ=0 ki=0 


kp-3=0 


^2^*^p-3'*^p-4* ‘ *'^0^ 


W 


T2k^^VN2 


0 


4 . z'-^. 


(3.5.51) 


where 
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To = 


ri-Ji 


2 ‘ 2 


= 2P-3j T + 2P'\_2 + ... + 2 j3 + J2 (3.5.52) 


and 


k(2) = 2P"^k..3 + 2P"^kp_4 + ... + 2ki + kp . 


In general, by letting 


9q(kp_q.] »l<p_q-2** * * »jq-2* * ** >jl 


1 

^ aq_1 (I<p.q9kp,q_i . ,k-| ,ko;jq_2,jq,3»- . . ,Ji >jo)Wo 
Kp-q-0 




j .|k(P)/N I , for q <>p 


W P 


(3.5.53) 


the general form of the Fourier coefficients^may^ bewritten as 


1 1 1 

H(r) = S Z S ®q(kp_q_i »kp_q.2» • • * >ki ,Icq; 


where 


ko=0 k-j=0 


>-q-l 


=0 


^q-1 ’^q-2’ ' ’ ' ’^1 


rqk(<)VNq 


(3.5.54) 


Nq - 2 


aii = 2P"^ 


^ Vi. - J .q.-l . 2P*'''^ , + zP'P-2j , + 

2 '^D-1 ^ ‘^P“2 ^ • 


and 


+ 2jq+i + jq 

k(q) = zP-^’-’k 1 + 2 P-P-Vq -2 + • • . + 2 ki + ko . 


(3.5.55) 


For the case p = q + 1, Eq. (3.5.54) may be written as 
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1 

H(r) = z 

ko=0 

where p ts the largest integer value of the time series with sample 
length N = 2^. Further reduction of Eq. (3,5.56) gives 

H(r) ~ 9p{ jp-] > jp-2» • * • *Jl »Jo^ (3.5.57) 

which is simply the transformation into a single data sample. 

3.6 Computing Techniques 

The computational procedures using the FFT method for power and 
cross spectra of multiple stationary time series Is based on the re- 
quirement that the total number of data points transformed to be 
integer powers of two. If the blocked time series do not have this 
required number of data points, zeros must be added to fill up the 
series to the required number before applying Fourier transformation. 

In case the calculated spectrum will be used to obtain the cor- 
relation function, then the time series should be filled up with zeros 
to obtain a total data points of to start with. This new approach 
of obtaining correlation functions by using two passes of Fourier 
transformation has been considered to be more economical in comparison 
with the direct multiply-add operations in obtaining the correlation 
functions. These correlation functions can further be used to obtain 
smoothed spectral estimates by applying the properly chosen lag 
windows whose characteristics are more generally discussed in time 
series analysis [51]. 


®p-i * ■ *■^1 




(3.5.56) 
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The general procedures used in computing power and cross spectral 
estimates via FFT are listed in the following two sections. 


3.6.1 Power Spectrum 

Reasonable power spectral estimates can be obtained for a time 
series of length n = 2^ by the following procedures: 

1. Either truncate the excessive data or add zeros so the number 

I 

of data points for blocked time series a^"^, i = 1,2 n, m = 1,2,..., 

M-1 , to be transformed is n = 2P, where p is an integer and M is the 
total number of blocks. 

2. Taper the blocked time series with a cosine taper data 
window as discussed in section 3.3.1 or with another appropriate taper- 
ing function as presented in section 3.3.2. 

3. Compute the finite Fourier transform of each blocked subseries 
by 

^ n-1 'in . 

H^^(f) = E at (t)W^'^ r = 0.1. 2,..., n-1 (3.6.1) 

^ k=0 


where W = e m = 1,2,..., M-1 is the block number. 

4. Compute the absolute squared value scaled appropriately to 
obtain the power spectral estimates by 



2At 

n 



r = 0,1,2,. ...n-1 (3.6.2) 


'' m ' m 

where subscript a in refers to blocked time series a^. for 1 = 

l,2,...,n and m - 1,2,...,M-1. 
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5. Adjust the power spectral estimates for a scale factor due 
to the cosine data window tapering by 

^ar “ Q Q 75 ^ar r=0*1*2*...,n-l . (3.6.3) 

6. Adjust the segment average for M-1 blocked spectral estimates 
each of length n by 

1 - m 

Gar = mT ^ r = 0.1,2.... ,n-l . (3.6.4) 

m=l 

7. Apply the frequency smoothing for segmently smoothed spectral 
estimates of length n = 8192. The spectral estimates are unique only 
up to the point where r « j + 1. This is due to the result of the 
application of the Fourier transform. At this point the Nyquist cut- 
off frequency occurs. The smoothing is performed for 4096 spectral 
estimates without including the zero frequency point. 

Looking at the printed-out spectral estimates of 4096 values, 
the energy is more concentrated in the lower frequency range. The 

proposed frequency smoothing technique as discussed in section 3.4.4 is 

therefore adopted. The frequency smoothed spectral estimates may be 
considered as representing the midpoint of the frequency interval from 
f|^ to fk+£-i, k ^ 1, where i is the chosen average value except when s, 
equals unity. The value of the frequency associated with different 
values of i and which corresponds to the spectral estimate in Eqs. 
(3.4.5) to (3.4.9) is calculated as follows 
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I = 1 


A = 4 


£ = 16 


i = 64 


fk = kAf 

fk = C^(k“7) + 2.5]Af 
fk = [I6{k-31) + 9.5]Af 


fk = [64(k-54) + 33.5]Af 

and z - 256 


k = 1,2.... .8 , 
k = 9, 10,. -.,38 , 
k = 39,40,... ,62 , 
k = 63,64,.. .,86 , 


(3.6.5) 

(3.6.6) 

(3.6.7) 

(3.6.8) 


fk = [256(k-79) + 129.5]Af k = 87, 88,..., 94 . (3.6.9) 

where Af = . The number of the spectral estimates was reduced to 

94 as a result of the application of the segment average and frequency 
smoothing technique in section 3.4.4. The frequency points represented 
In Equations (3.6.5) to (3.6.9) are corresponding to these smoothed 
spectral estimates. 


3.6.2 Cross Spectrum 

The computations for the cross spectral estimates between two 
discrete time series a^ and b-j , i = 0,1 ,2,. . , ,n-l seem to be more in- 
volved due to the requirement of obtaining the coherence function and 
phase angles between the two time series. The general procedures for 
calculating cross spectral estimates via FFT can be written as follows 

1. Either truncate the excessive data or add enough zeros in 
the two discrete time series and b^-^, i = 0,1 ,2,. . . ,n-l so that 
n = 2P as discussed in step 1 in section 3.6.1. 
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2. Taper both blocked time series and i = 0,1 ,2, . . . ,n-l 
with the cosine taper data window or other appropriate data windows 

as presented in section 3.3.2. 

3. Let be the real part and b-j*^ be the imaginary part of a 
newly created complex time series c^'*' = a*"" + jb^-'^. i = 0,1 .2,. . . ,n-l 
and j = /T . 


4. Compute the finite Fourier transform of each blocked complex 
time series by 


Cr'^(f) = 2 (t) W'"'^ r = 0,1,2 n-1 (3.6.10) 

- 2tt ^”0 

where W = e and m = 1,2,...,M-1 is the block number. 

5. Obtain the respective spectral values of both a^"^ and 
'm 

^i > i • 0,1 ,2,. . . ,n-l by using the properties of the Fourier trans- 
formation as discussed in section 3.5.1 as follows 

. C'r'"(f) + CrW (3.6.11) 

‘'n 2 


and 


•V m, c 


Af) 


2j 


Cn-r(f) 


(3.6.12) 


•^m ~m* 

where Cn.y,(f) represents complex conjugate of Cn_y.(f) and j = /T. 

6. Compute the cross spectral estimate r = 0,1 ,2,. , . ,n-l for 

block m by 


- m 2At 
•3abr = n 


-m* 

pin 

tr-j 


m 


(f) Cri(f) 


(3.6.13) 


7. Adjust the cross spectral estimates by a scale factor due to 
cosine data window tapering in each block 
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A m A in 1 n 

Gabr(f) = Gabr(^^) r = 0,1 ,2, . . .^ + 1 (3.6.14) 

8 , Apply the segment average smoothing for M-1 blocks of spectral 

estimates each of length n by 

A 1 M-1 . m 

^abr^^) “ M-T ^ *^abr r = 0,l,2 j+1 (3.6.15) 

m=l 

9. Pick up the real and Imaginary part of 6 ^ 5 ^, by 

Gabr(f) = Coabr(f) - 3 Qabr^^') >• = 0.1 .2, . . .| + 1 (3.6.16) 

where Co^brff) called the cospectral density function and Qabrl"*") 
is called the quadrature spectral density function. 

10, In order to apply the proposed frequency smoothing technique 
for cross spectral density estimate, two different ways can be taken 
as the following 

a. Smooth the cross spectral estimate Ggbr(^) ob- 

tain the Co- and quadrature-spectral density function. 

b. Smooth the real and imaginary parts of Gaby.{f) individually 
and then obtain the final smoothed cross spectral estimates. 

These two prcoedures can give different results, since the 
linear smoothing operation, the non-linear operations of the square 
root of squared value, and division are not corrinutative. In meteoro- 
logical applications, it is the co- and quadrature-spectral density fu- 
nctions that are of interest instead of the positively valued cross 
spectral density function. Thus, the procedure of smoothing the real 
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and imaginary parts of the cross-spectral estimates separately is 
used in this case. The final smoothed cross spectral density esti- 
mates are obtained by 

Gabr(f) = ■^COag^(f) + r = 1 .2, . . . ,94 . (3.6.17) 


n. The phase angle contained in the smoothed cross spectral 
estimate Gabr('^) calculated by 


®abr “ 


Qabr(f) 

_Coabr(^). 


r = 1.2 94 . 


(3.6.18) 


12, The smoothed squared coherence estimate is given by 



(3.6.19) 


where Ga^(f) and Gbr(f) are calculated by Eq. (3.6.4). The frequency 
smoothing technique of Gay.(f) and Gb^(f) was presented in the previous 


section. 



CHAPTER IV 


STATISTICAL ERRORS 

Errors in calculations of statistical quantities of digitized 
time series are quite uncertain because of the large amount of data 
collected* the underlying probabilistic nature of the data and the 
method in deriving the desired statistical parameters. 

The random nature of the data makes it almost impossible to know 
the deterministic characteristics of a physical phenomenon. It is 
only possible to know the average level and to obtain some estimate 
of the reliability or accuracy of this average level. 

In representing an estimate of a statistical parameter* no indi- 
cation of the reliability of the estimate is found from the simple 
calculations of the estimated value. Given the size of the sample does 
not provide the means to interpret the accuracy of the estimate as 
a function of the sample size. A more direct indication of the 
accuracy is desirable. 

The standard error of any estimate is used to indicate the 
reliability or precision of the calculation of this estimate. What 
one really wants is a range of values within which the estimates 
should fall. Therefore, the notion of confidence interval for a 
parameter is commonly used to serve this purpose. 
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This chapter will investigate the validity of estimating the 
blocked mean and blocked variances in the blocked time series in 
relation to the sample mean and sample variances of the total sample. 
The confidence interval of spectral density estimates as a result of 
using the FFT technique will be determined in this chapter. 

4.1 Mean Value Estimate 

The block mean value for n data points of the i^^ time series in 
the mth block may be calculated by 

Ai™ = i Z Ai(j) i =1,2 4 (4.1.1) 

where bar denotes the mean value and i indicates the number of the 
time series. 

The sample mean of the time series is simply the arithmetic 
average of the block mean values as obtained fay 
_ 1 M _ni 

Ai = ]([ S Ai i = 1,2.... .4 (4.1.2) 

m=l 


where M is the total number of blocks in the time series. 

In case the time series has not been blocked, then the sample 
mean is calculated by 


_* 1 ^ 

Ai “ N ^ A<j ( j ) 
« j=l 


i 


1.2 i 


(4.1.3) 


where N is the total number of samples in the time series (i.e. N = 
nM). 
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Combining Eqs. 


Mn 


m=l j=l 


(4.1.1) and (4.1.2), one obtains 
At(j) 


1_ 

Mn 


n 

S Ai(j) 

j=l 


2n 

+ I A,-(j) 
j=n+l 


Mn 

+ s A^-(j) 
j=(M-l)n+l 


1 Mn 

" Mn ^i^^^ 
J ^ 


i = 1 , 2 ,..., 4 


(4.1.4) 


Comparing Eqs. (4.1.3) and (4.1.4), it is seen that the sample 
mean of the sample can be calculated from the average of the block 
mean values in each subseries. The only possible difference between 
the calculated sample mean value of Eq. (4.1.2) and (4.1.3) is the 
summation of large numbers by floating-point representation in the 
computer without using double precision. In fact, the sample mean 
values calculated from Equation 4.1.2 should be more accurate than 
using Equation 4.1.3 directly. 


4.2 Variance and Covariance Estimates 

The sample variance and covariance are calculated from the fil- 
tered series with exactly zero mean by Eq. (2.4.2) as follows 

1 M-1 

T 

m=l 


®i^j = 


i-1 ®i®j 


-m 


i.j = 1 »2 4 


(4.2.1) 


where block variances 


^i^j 


-m 


n 

I 

k=l 


and covariances are calculated from Eq. 
a^(k)aj(k) i , j = 1,2,..., 4 


(2.4.1) 

(4.2.2) 
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where a-j, i = represents filtered time series. Equal sub- 

scripts i “ j in Eq. (4.2.1) denote the variances while unequal sub- 
scripts i t j represent the covariances. 

In case the filtered subseries has a near zero mean but not 

exactly zero mean, then the values of block variances and covariances 

calculated by Eq. (4.2.2) will be overestimated. The variance and 

covariance of the total sample with near zero mean values for the sub 

series are calculated as follows 

* 1 M-1 1 n ^ ■ 

®i®j M-1 n (af(k) - a,* ) (aj(k) - aj ) 

1 M-l r, n ^ n 

n ®iO<)ajOO 





M-1 

s 

m==l 





M-1 


M-1 

I 

m=l 


— m — 
^1 


(4.2.3) 


where 

— m 1 ^ . 

a,- = n X a^(k) 1 = 1,2 4 . (4.2.4) 

k“ 1 

IP 

Consequently, a^ is the mean value of the subseries after the moving 
average filter has been applied to the total sample. The reason we 
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chose the number of blocks to be M-1 here is because the data samples 
were filtered and the ends were lost. Now, the sample variances and 
covariances based on the total filtered sample can be calculated as 
follows 

1 N __ 

= N 2 : (a^(k) - a^-) (aj(k) - aj) , (4.2.5) 

k“ 1 

where aY is defined by 

— 1 N 

a. = M E a,.(k) 1 = 1,2,... ,4 (4.2.6) 

^ k=l ’ 

and N is the total number of data points in the filtered sample. In 
this case, N is equal to n(M-l) data samples. 

Expanding the terms in the right. hand side of Eq. (4.2.5) and 
combining, one obtains 

= a^-aj - a.j aj i,j = 1,2 4 (4.2.7) 

where the Equations 4.2.1 and 4.2.2 have been used. 

The error E(m) in variances and covariances as a result of the 
subdivision of the original time series in blocks may be obtained by 
combining the Equations 4.2.3 and 4.2.7 as follows 

i t 

E(m) = a^.a. - a.aj 

1 — tn — ra /» n o\ 

= M-1 I a,- aj - a.j aj (4.2.8) 


or 
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" M-1 ^ ^ • (4.2.9) 

m=l ^ 

The value of this error E(m) can be discussed by the following 
two limiting cases: 

1. In the limit when m = 1* the total sample consists of one 
block only, the error is zero. 

2 . In the limit when m approaches infinity and the sample con- 
tains a finite number of data points, then each block contains only 
one data point. In this case, the block mean value is simply the 
value of the data point itself in the block. It is obvious that the 
error E(m) is again zero. 

In case when the number of blocks is chosen to be a finite value 
as in most practical situations, the error E(m) can be significant 
since it depends on the time of the subseries. The time of the sub- 
series will affect the shape of the moving average filter as can be 
seen from Figure 5. The effective cut-off frequency in the filter with 
different time Interval is different. One should choose the block 
length of the subseries to be comparable to the period of the minimum 
frequencies of interest. The storage limitation in the computer for 
numerical calculations imposed another restriction in the selection 
of the block length in the subseries. 

4.3 Spectral Density Estimate 

For a stationary Gaussian random process {a(t)>, the power spectrum 
for a sample record a|^(t), k=0,l,.,n-l of finite length n is defined by 
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Ga(f) = 2 iim [|H(f,n)|2] (4.3.1) 

where 

H(f,n) = At z a^Ct) . (4.3.2) 

k=0 


The estimate of power spectrum can be obtained by omitting the 
limiting and expectation operations in the Eq. (4.3.1) to yield the 
following 


= -h. 


(4.3.3) 


with the narrowest possible resolution Bg = . 

For the discrete frequency values, the Fourier components are 
defined by 


„ H(f.n) 


n-1 

2 \{t) 

k=0 


^-j2irfkAt 


(4.3.4) 


and the spectral estimate is given by 

Ga(f) = ^r 


(4.3.5) 


The estimate of the spectral 
Gg(f) and is unbiased if 

E[Ga(f)] = Ga(f) 

The mean square error of the 
defined by 


density function Ga(f) is denoted by 

(4.3.6) 

power spectral density estimate is 


m.s.e. = E[(Gg(f) - Ga(f))^] . 


(4.3.7) 



In general, the power spectral density estimate is a function of 
the number of data samples n. In order that &a{f) is the consistent 
estimate of the power spectral density function Ga(f), it is required 
that 

£im E[(Ga(f) - = 0 . (4.3.8) 
n^ 


By assuming the random process {a.|}, i = l,2,...,n to be un- 
correlated white noise. Otnes and Enochson [67] proved that the power 
spectral density estimate Gia(i') an unbiased but not a consistent 
estimate of the power spectral density function 6a(f). 

4.3.1 Chi-Square Distribution 

The Fourier components H(f,n) computed by Eq. (4.3.2) are complex 
with real and imaginary parts, Hp(f,n) and Hj(f,n) which are uncor- 
related random variables with zero n^an and unit variance [6]. Both 
HR{f,n) and H (f,n) will be Gaussian random variables if the data 
sample a^•(t), i » 0,1 ,2, . . , ,n-l is Gaussian as a result of the linear 
operation of the Fourier transformation. It is seen that the quantity 

= HR2(f,n) + Hj2(f^n) (4.3.9) 

is the sum of the squares of two independent Gaussian variables. 

From the definition of Chi-square distribution, each frequency com- 
ponent of the power spectral density Gg(f) will have a sampling dis- 
tribution given by 
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2 

where X 2 the Chi-square variable with two degrees of freedom. The 
number of degrees of freedom, y, represents the number of independent 
or "free" squares summing in the expression (4.3.9). The mean and 
variance of the chi-square variable are y and 2 , respectively. The 
normalized mean square error can be obtained by 


e 2 - - 6a(f))^] 2 

[Ga(f)]2 ■ ^ 


(4.3.11) 


where use has been made of Equation 4.3.6. 

For Y = 2, the normalized standard error is given by 

±1 (4.3.12) 

which means that the standard deviation of the pSD estimate is as 
great as the quantity being estimated. To reduce the error of the pSD 
estimate as calculated by Eq. (4.3.5), smoothing the estimates is 
necessary. By smoothing either segmental ly or frequency smoothing 
as discussed in chapter three, the number of degrees of freedom y can 
be increased. It is seen, from Eq. (4.3.12) that the normalized 
standard error can be reduced if y is increased. 

The (1 - a) confidence interval for the power spectral density 
function Ga(f) around the frequency f based upon an estimate Ga(f) 
measured with a resolution bandwidth Bg and a record length n is 
given by 


VGa(f) 


2 

X yOL 


£ Ga(f) 


YSa(f) 




< 


- 1 - a . 


(4.3.12) 



89 


and Y = 2Bgn . (4.3.13) 

The term, 1 - a, is a fixed confidence level which commonly is taken 
to be 0.80, 0.90 or 0.95. The true value 6g(f) lies between the two 
values in the bracket of Eq. (4.3.12). The Chi-square distribution 
Ts tabulated in Table III and defined by 

= [b such that f P(x^)dx? = a] (4.3.14) 

v;ot r 'y 

For degrees of freedom where y > 30, the following expression 
may be used to obtain the distributions. 

^1-a * ? (4.3.15) 

* 

where T^ is the corresponding percentile of the standard normal dis- 
tribution (Table I). 

4.3.2 Numerical Example 

An example will illustrate the application of Equations 4.3.12 
and 4.3.15. 

The total number of data points used for the estimation of the 
spectral density for 43 blocks each of which contains 8192 data points 
is 

N = 8192 X 43 = 352.256 . 

The effective bandwidth for a sampling rate of 200 samples per second 
is 
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Be = fj^ = 0.000567 HZ . 

The number of the degrees of freedom is 
Y - 2BgT == 2 . 

The 95% confidence interval for two degrees of freedom in Chi-square 
distribution is found from Table 3 to be 


x|; .975 = -051 
4 -, .025 = 


(4.3.16) 


A 

Ga(f) is the spectral density estimate, then the 95% confidence 
interval of the true spectral density function Ga(f) is given by 

7^ Ga(f) < Ga(f) < q7o 51 ^a(^)] 
or 

[0.271 Ga(f) < Ga(f) <39.215 Ga(f)] . (4.3.17) 

It is seen from Eq. (4.3.17) that the true spectrum lies in a 
wide range of values and that little confidence can be placed on the 
estimate. The range of the confidence interval can be reduced by in- 
creasing the number of degrees of freedom. The number of degrees of 
freedom can be increased greatly when the segment average technique is 
applied. By averaging over 43 blocks of spectral density estimates, 
the resulted smoothed spectral density estimate has the number of 
degrees of freedom as follows 
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Y = 43 X 2 = 86 

The 95% confidence bands for 86 degrees of freedom in Ch1-square dis- 
tribution is formed by using Equation 4.3.15 as follows 


4e; .975 = 61.68 

><86; .025 = H3.2 


(4.3.18) 


and the true spectrum lies in the interval 

[.76 Ga(f) < Ga(f) <1.39 G^tf)] (4.3.19) 

It is seen that the range of values in between which the true spec- 
trum lies has been reduced as a result of applying the segment smoothing 
technique in the spectral density estimates. The application of 
frequency smoothing on the spectral density estimates may give a much 
better representation in the range of estimated values. 

The analysis shown above was based upon the assumption that the 
data samples were both Gaussian and white. In general, the data are 
colored (i.e. correlated) in some manner. This has the effect of 
reducing the number of the degrees of freedom in the Chi-square dis- 
tribution. The standard practice is to use the white noise results 
as a guideline 1n the spectral density estimates [67]. 



CHAPTER V 


DISCUSSION OF THE RESULTS 

The data used for the statistical analysis in this dissertation 
were measured with either the Model 1080D Total Vector Anemometer 
probe {triple split films) or the Model 1296L dual split film 
probe both manufactured by Thermo-Systems, Inc. The operation and 
the analysis of the data from these probes are discussed in detail 
in references [95] and [108] respectively. In conjunction with the 
two types of TSI probes a set of Gill propeller anemometers was used. 
These propeller anemometers were mounted in such a way that one was 
parallel to the TSI probe and the other perpendicular to the first 
one and both in a horizontal plane and located adjacent to the TSI 
probe. A detailed discussion of the operation of the Gill anemometers 
and data analysis can be found in Appendix A. The Gill anemometers 
were used mainly for comparison of results with those obtained from 
both types of TSI probes. The anemometers were located on the top of 
the air exchanger of the low speed wind tunnel at Virginia Polytechnic 
Institute and State University. This was the best location available 
in the neighborhood of this wind tunnel in which the anemometers were 
calibrated. The connecting cable between the probe and the anemo- 
meter was 350 feet long, and since it was not feasible to move the 
trailer in which the data acquisition system was located, the above 


92 


93 


described location for the probes was the best available within 
a 350-foot radius from the trailer. 

It turned out that the winds measured in this location as mentioned 
above were mainly from a north-west direction. However, due to the 
presence of upstream buildings and due to the location of the anemo- 
meters on top of the wind tunnel a flow was measured with an appreciable 
vertical component and with relatively large fluctuations in magnitude 
and lateral direction. As a result of the type of data measured, a 
great deal of effort was used in the proper analysis of the data with 
respect to stationarity, filtering and smoothing of the calculated 
statistical quantities. Mean values, variances, covariances, power 
spectral estimates, cross spectral estimates, coherence functions and 
phase angles of the three turbulent wind components in the mean wind 
oriented coordinate system were calculated. 

During the period in time that these data were taken, considerable 
troubles were still encountered with the data acquisition system 
(see reference [95] for a detailed description of the data acquisition 
system). Specifically, the PDP-11 /20 mini -computer had an intermittent 
problem which was extremely difficult to pin down. Also, the data 
acquisition as well as the consequent digitizing of the data was 
affected by the simultaneous operation of the wind tunnel. It was 
found that fluctuations in line voltage as a result of the starting of 
the wind tunnel often erased the recorded time of day or the recorded 
voltages of the anemometers. 
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Due to the fact that two different anemometers (TSI and Gill) 
were used simultaneously, the data from the Gill had to be recorded 
on an F.M. tape recorder first and consequently played back at a 
later time in order to be digitized. The data received from either 
TSI - anemometer was fed directly into the multiplexer and directly 
stored on the digital tape. It often occurred that bad data or no 
data were recorded on the digital tape or that tape marks disappeared 
and as a result the data were useless. Since this period, the entire 
data acquisition and data handling system has been moved to NASA - 
Wallops Station for data gathering from the meteorological tower at 
Wallops Island. The Wallops Island data will be analyzed at a later 
date, only data taken at the Virginia Polytechnic Institute and State 
University location will be discussed in this dissertation. The data 
acquisition and data handling system while in operation at Wallops 
Island did not experience the breakdown and problems as were en- 
countered when in operation at Virginia Polytechnic Institute and 
State University. Consequently, the data discussed in this report 
are somewhat sketchy, but adequate to indicate that very good 
results can be obtained with the system. 

The accuracy and the efficiency of the computer program which 
was developed in order to calculate the statistical quantities of a 
digitized time series, was tested against a simulated time series. The 
theory and the processing of digital simulation of random processes 
is described in detail by Sinha [87]. If the power-spectral density 
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function is known, a sample of the corresponding time series can be 
developed by inverse fast Fourier transformation. 

In order to simulate the streamwise turbulence component in 
the atmosphere, use was made of the semi -empirical von Karman 
spectrum given by the following expression in dimensionless form 


fG(f) ^ 


4U 

U 


+ 70.78 


A 20% intensity of the turbulence was assumed so that = 0.20 IT. 
Here, G(f) is the power spectral density function of the streamwise 
turbulence component so that 

‘OD 

6(f)df = ? , 

*^0 


and L is the longitudinal integral scale which varies with height but 

I 

which was chosen to be 360 feet. 

Based on this information a value for the simulated time series 
was calculated at intervals of .05 seconds. In this manner (2)^^ 
data points were generated which represents the digitized time series 
sampled at a rate of 20 samples per second and of 27.31 minutes 
duration (Figure 12). This simulated, digitized time series was then 
used as the input to the computer program which was developed for the 
statistical analysis of time series of long duration. The calculated 
power spectral density function was then compared with the original 
semi-empirical von Karman spectrum from which the time series was 
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developed. As can be seen from Figure 13, the similarity is very 
good, especially at higher frequencies. At lower frequencies, the 
estimates of the power spectral density function seem to deviate a 
great deal from the expected values which is due to the limited 
sample length of the simulated time series. 

Because of the way the fast Fourier transformation works, only 
a few points at the lower frequency range are calculated and as a 
result one can only resort to the so-called segment smoothing of 
the spectral estimates in this range. At the high frequency range 
the density of calculated spectral estimates is much higher and as 
a result the so-called combined smoothing process can be applied 
with the result that the estimates show much less scattering. 

The statistical quantities are calculated from the data 
measured by three different wind measuring sensors, namely, Gill 
anemometers, TSI Model 1080-D total vector anemometer and the TSI 
Model 1296, dual split-film probe for four separate runs. Mean 
values, variances and the covariances are all listed and compared in 
Tables IV through VII. The smoothed spectral density estimates are 
plotted versus frequencies in Figures 14 through 25. 

In Table IV, the block means, samples means and the number of 
reverse arrangements of the block means used for statistical test are 
listed for each velocity component in the sensor oriented coordinate 
system. Mean wind velocity, variances and covariances of the ve- 
locity components and temperature in the mean wind coordinate system 
are also tabulated. The calculated numbers of reverse arrangements of 
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the block means of the velocity components show that the hypothesis of 
nonstationarity as far as the trend in the data is concerned is accept- 
able at the 10 percent level of significance for the A and B components 
of the veloicty. The reverse arrangements of the block means of the C 
component shows that some trend is present. However, for this run the 
results of the trend test are in general acceptable. If the number of 
reverse arrangements would be significantly different from the expected 
numbers given in Table II, the hypothesis of acceptable nonstationarity 
would be rejected and the data would be rejected for further analysis. 

The results of run 10 show that an appreciable vertical mean ve- 
locity component is present and the turbulence intensities are between 
20 and 25 percent for all three components. 

The smoothed power spectral density estimates of each velocity com- 
ponent of run 10 are plotted against frequency in Figures 14 through 16. 

The area under the curves correspond closely to the respective variances 
and at the high frequency range the spectrum functions vary as expected 
as the frequency to the -5/3 power. Due to the limitation of core storage 
In the IBM 360/155 digital computer, the number of data points selected 
in the subseries restricts the computation of the spectrum functions at 
frequencies above 0.0244 hertz (period of 40.96 seconds). Due to the 
chosen sample rate of 200 samples per second the maximum frequency at which 
the spectrum function can be analyzed is 100 hertz. 

The cospectrum between the longitudinal and vertical velocity com- 
ponents is plotted in Figure 17, and the shape of the spectrum agrees 
fairly well with that which was observed at Brookhaven, Long Island and 
and reported by Panfsky [72]. 
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In Table V, the results calculated from the data of run 11 as 
measured by the TSI #1192 total vector anemometer and the Gill ane- 
mometers are compared. In order that the results calculated from the 
data as measured by these two different sensors are comparable, the 
block means calculated from the velocity components measured by TSI 
#1192 probe are transformed into the probe oriented coordinate system. 

The number of reverse arrangements of the block means for measurements 
from the TSI #1192 probe are calculated based on the velocity compo- 
nents in the sensor oriented coordinate system. The number of reverse 
arrangements of the block means for Gill anemometers are calculated 
based on the data obtained In the probe oriented coordinate system. 

The number of reverse arrangements for the block means of the 
lateral velocity component is very high and after examination of the 
data one can see that this is due to a gradual change in wind direction 
or due to the gradual change from positive to negative lateral velocity 
components. The longitudinal conponent of the velocity does nbt show 
any trend and is not very much affected by the gradual change in wind 
direction. Comparison of the block means and the sample mean shows ihat 
the longitudinal component measured by the Gill anemometer is consistently 
higher than that measured by the TSI probe. The block means and the 
sample mean of the lateral velocity component and the angle the mean 
wind makes with the direction of the instruments compare very well. 

Due to the limited response characteristics of the Gill anemometers, 

the turbulence quantities obtained from the Gill anemometer are consistently 

lower than those obtained from the TSI probe. Due to the presence of some 
high frequency noise during the data runs, the smoothed spectral estimates 
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have folded in the high-frequency range and as a result the graphs for 
the smoothed spectral estimates for run 11 are not shown. 

In Tables VI and VII, the results calculated from the data of 
both runs 12 and 14 as measured by the dual split-film anemometer and 
Gill anemometers are tabulated and compared. Both these anemometers 
measure the velocity components In the same probe oriented coordinate 
system. Since the shafts of the propellers of the Gill anemometers 
are parallel to each of the two perpendicular coordinate axes in the 
horizontal plane, only the block means of the longitudinal and lateral 
velocity components can be compared in Tables VI and VII. The block 
means calculated from the vertical velocity components as measured by 
the dual split-film probe are tabulated only for completeness not for 
comparison. 

The number of reverse arrangements for the block means of the 
longitudinal and lateral velocity components for both run 12 and 14 show 
that the data are free of abnormal trends. Again the Gill anemometer 
overestimates the block means of the longitudinal velocity component 
by as much as 15 percent. The variances \? and are estimated lower 
by the Gill anemometers as before. For this particular type of flow the 
values for the covariance uw are very low and compare reasonably well. 
The horizontal angles of attack between the mean wind and the instruments 
compare quite well for both runs. 

In both Tables VI and VII, comparisons are also made for the values 
of the sample means, the variances, and the covariances of the velocity 
components in the mean wind coordinate system for the cases with either 
the near-zero block-mean values removed or not removed. These near-zero 
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block means are the mean values remaininci In the data samples of each 
block after having applied the moving-average and differencing high- 
pass filter to the origional time series as discussed in chapter two. 

The small near-zero mean values do not have any effect upon the values 
of the mean velocity components in either the probe oriented coordinate 
system or the mean wind coordinate system since the mean velocity com- 
ponents are calculated from the unfiltered original time series. The 
effects of these near- zero means upon the estimated values of the 
variances and covariances in the mean wind direction are not very 
significant. They do have significant effects upon the estimated values 
of the spectral density function at low frequencies. 

Figures 18 through 20 show the plots of the power spectrum of 

the longitudinal, lateral, and vertical velocity components measured by 
TSI #122 in run 12, respectively. The power spectrum of the longitudinal 
and lateral velocity components measured by the Gill anemometers in the 
same run 12 are plotted in Figures 21 and 22. It is quite evident when 
like spectra form the TSI probe and the Gill anemometers are superimposed, 
the difference at frequencies higher than one hertz are considerable! 
expecially when one realizes that the ordinate of these spectrum functions 
has a logarithmic scale. For both spectrum functions of the longitudinal 
as well as the lateral turbulence components the Gill anemometers show a 
too rapid decrease with frequency in the range from 0.5 to 100 hertz. In 
the low frequency range the comparison is quite good. 

The power spectrum of the longitudinal velocity component measured 
by the TSI #122 probe in run 14 is plotted in Figure 23. In order to 
obtain some spectral properties in the lower frequency range below 0.024 
hertz, two new time series were generated from the data samples in run 14. 
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The first one (series A) is generated by taking a 10-point non-over- 
lapping average from the original time series. The second one (series 
B) is obtained by containing every 10th point from the original time 
series. The two newly generated time series have a reduced number of 
data points to one-tenth of the original number. The sampling rate 
is also decreased from 200 sanples per second to 20 samples per second. 

The spectral density estimates for the newly generated time series are 
calculated by following the similar procedures used in computing the 
estimates from the original time series except the segment average is 
taken over the corresponding estimates of the four blocked subseries. 

The power spectra calculated from the time series A and B are 
plotted in the Figures 25 and 24 respectively and can be compared with 
the original power spectrum of run 14 plotted in Figure 23. The 
power spectra of the newly created time series show a greater deal of 
scatter in the low and intermediate frequency range. The spectrum of 
the time series B shows some frequency folding in the high-frequency 
range since the data were not filtered at 10 hertz. The spectrum of 
the time series does not show this folding since the averaging 
procedure act as a low-pass filter at approximately 10 hertz. At the 
low-frequency range the scatter of the spectrum data for both time 
series A and B is more severe than for the spectrum data of the original 
time series. The limited segment averaging results in a reduction of 
the number of degrees of freedom which will give a broad confidence in- 
terval In the spectral density estimates. The data in each block sampled 
at a rate of 20 samples per second allows calculation of spectral estimates 
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down to a frequency of 0.00244 hertz (period of 409.6 seconds) since 

1 3 

the number of data per block is limited to 2 or 8192. However, only 
four blocks can be obtained from time series A and B,and in the low- 
frequency range where only segment averaging can be applied the spectral 
data can only be averaged over four values as compared to 44 values 
in the original time series. In order to get more reliable spectrum 
estimates at frequencies lower than 0.0244 hertz, the total sample time 
should be increased ten times if spectral information down to 0.00244 
hertz is required. However, with time series A and B we can get some 
idea how the spectrum varies at frequencies lower than 0.0244 hertz. 


CHAPTER VI 


CONCLUSIONS 

The statistical analysis of a discrete time series with a large 
number of data points representing a random process can be achieved 
successfully by the procedures as outlined in this report. It is 
necessary that these time series are subdivided into a certain number 
of sample records each containing an equal number of data points. The 
number of data points in each sample record depends largely on the 
storage capacity of the available digital computer. The total number of 
data points in each sample record should be chosen such that it is an 
integer power of two in order to satisfy the requirements for the fast 
Fourier Transformation. The statistical quantities calculated from 
each sample record or data block can be used to determine the degree of 
stationarity of the total sample by application of the nonparametric 
statistical test. The existence of any trends in the time series can be 
removed successfully by using the moving-average and differencing high- 
pass filter. This type of statistical analysis was used to obtain 
statistical information from long time series representing low-level 
atmospheric winds and temperature. These quantities were measured with 

fast response split-film anemometers developed by Thermo-Systems Inc. 

1 

A set of propeller-type Gill anemometers was used simultaneously to 
measure wind velocities in the horizontal plane in order to compare their 
results with the results from the TSI probes. 

The digitized data obtained from the TSI probes are stored on 
digital tape as voltages. Seven voltages are necessary to obtain three 
velocity components and temperature if the triple split-film probe is 
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used, and five voltages are necessary to obtain the same information 
when the dual split-film probe is used. The following steps are re- 
quired to obtain the statistical information for each data set. 

A. The first step is to calculate for each sample point the three 
velocity components and the terrperature. This information is obtained 
through calculations as outlined in reference 95, and 1t is consequently 
stored on another digital tape In blocks of 418 sample points each. The 
velocity components are calculated in the so-called sensor-oriented co- 
ordinate system. 

B. In the second step, the block size is changed from 418 sample 
points to 2 ''^ = 8192 sample or data points. First the meansof each ve- 
locity component and the magnitude of the mean velocity for each block 
are calculated. Also the mean temperature and the average horizontal 
angle the mean velocity for each block makes with the probe axis is 
obtained. In addition, the block variances and covariances for the three 
velocity components and the temperature are calculated. The number of 
blocks depends on the recording length of the data, but for a run of 
about one half hour the number of blocks is 44. In addition to the block 
means, the sample means of the velocity components, the temperature as 
well as the horizontal angle between the sample mean-wind and the axis 

of the probe are calculated. Also the number of reverse arrangements of 
the block means of the velocity components and the temperature and those 
of the block standard deviations of velocity components and temperature 
are calculated. The latter calculations are made in order to check for 
nonstationarities such as time-varying mean values or time-varying standard 
deviations or a combination of these two. At this point further analysis 
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of the data would be stopped if the number of reverse arrangements would 
be significantly different from the expected numbers as given in Table II. 

C. The third step of the data analysis consists of the removal of 
the mean and low frequency components of the four time series. This is 
accomplished with the use of the moving-average and differencing high- 
pass filter. After this filter has been applied to the data four new 
time series with new sample mean and without low-frequency components are 
established. The cutoff frequency for this type of filter with a filtering 
interval of 40.96 seconds (the length of one block) is approximately 

.0108 hertz. 

D. In the fourth step the sample variances and covariances of the 
four filtered time series are calculated and consequently transformed 
into the mean-wind coordinate system. 

E. In this step the data representing the filtered velocity com- 
ponents in the sensor oriented coordinate system are transformed into 
components of the mean-wind coordinate system. As a result four time 
series with zero mean are created representing the fluctuating temperature 
and, the fluctuating velocity components in the mean-wind coordinate 
system. 

F. In this step the power spectral densities of the four time 
series obtained in step E are calculated using the newly developed fast 
Fourier transform method with the no-bit reversal procedure and also using 
the appropriate combined smoothing techniques. 

G. In the last step the coincident spectral density function and 
the quadrature spectral density function of two different time series 
are calculated using the newly developed fast Fourier transform method 

and in addition using both the segment as well as frequency averaging method 
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In order to check the accuracy of the calculated power spectral 
density extimates of any time series, spectral densities were calcul- 
ated from a simulated time series which was generated from a known 
spectrum function. The calculated spectral densities correspond quite 
well with the spectrum function from which the simulated time series 
was generated. This indicates that with the newly developed computer 
program using the fast Fourier transform with a no-bit reversal procedure 
and with proper smoothing procedures, accurate spectral information can 
be obtained in the frequency range between 0.0244 and 100 hertz. Spectral 
density estimates of a lower degree of accuracy for frequencies less 
than 0.0244 hertz can be obtained by creating a new time series by taking 
as an exanple a 10-point non-overlapping average of the original time 
series. Mean values, variances as well as covariances of the two hori- 
zontal wind components measured with the TSI probe were compared with the 
same quantities measured simultaneously with the Gill propellor anemometers. 
The discrepancies in these quantities can be attributed to the varying 
and limiting response characteristics of the Gill propeller anemometers. 

As a result of carefully carried out calibration procedures and 
the application of the newly developed computer program, accurate statis- 
tical estimates of long time series describing the fluctuating wind com- 
ponents can be obtained with either the dual or the triple split film TSI 
probes. 
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Figure 2. Effectively Removable Trend Types for the Mean 
Square Value. 
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Figure 6. The Effect of Slow Sampling Rate of a High Frequency Wave. 




127 


TRIANGULAR SPECTRAL WINDOW 




Figure 9. Filter Shape Before (top) and After 
(bottom) Segment Smoothing. 
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INITIAL TRIANGULAR SPECTRAL WINDOW 



Figure 10. Filter Shape Before (top) and After (bottom) 
Frequency Smoothing (K Is the number of 
smoothing points). 
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Figure 11. Overall Filter Spacing of the Proposed Frequency Smoothing Technique 
(m = 8192, At = 0,005 second). 
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Figure 13. Power Spectrum of the Simulated Time Series 
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Figure 14. Power Spectrum of the Longitudinal Velocity Component. 
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Figure 15. Power Spectrum of the Lateral Velocity Component. 
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Figure 16. Power Spectrum of the Vertical Velocity Component. 
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Figure 17. Cospectrum Between Longitudinal and Vertical Velocity 
Components. 
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Figure 20. Power Spectrum of the Vertical Velocity Component. 
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Figure 21. Power Spectrum of the Longitudinal Velocity Component. 
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Figure 22. Power Spectrum of the Lateral Velocity Component. 
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Figure 23. Power Spectrum of the Longitudinal Velocity Component. 
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Figure 24. Power Spectrum of the Longitudinal Velocity Component. 
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TABLE I. PERCENTILES 0 
DISTRIBUTION. 


P[T < T^] 

* 

Ta 

.001 

-3.09 

,005 

-2.58 

.010 

-2.33 

.020 

-2.05 

.025 

-1.97 

.030 

-1.88 

.040 

-1.75 

.050 

-1.645 

.100 

-1.28 

.150 

-1.04 

.200 

-0.84 

.300 

-0.52 

.400 

-0.25 

.500 

0 
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THE STANDARD NORMAL 


T < T’] 

T 

‘ct 

.600 

0.25 

.700 

0.52 

.800 

0.84 

.850 

1.04 

.900 

1.28 

.950 

1.645 

.960 

1.75 

.970 

■1.88 

.975 

1.97 

.980 

2.05 

.990 

2.33 

.995 

2.58 

.999 

3.09 
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TABLE II. PERCENTAGE POINTS OF REVERSE ARRANGEMENT DISTRIBUTION.* 


£ « 



0.99 

0.975 

0.95 

0.05 

0.025 

0.01 

10 

9 

11 

13 

31 

33 

35 

n 

12 

14 

16 

38 

40 

42 

12 

16 

18 

21 

44 

47 

49 

13 

19 

22 

25 

52 

55 

58 

14 

24 

27 

30 

60 

63 

66 

15 

28 

32 

35 

69 

72 

76 

16 

34 

38 

41 

78 

81 

85 

17 

39 

44 

48 

87 

91 

96 

18 

45 

50 

54 

98 

102 

107 

19 

52 

57 

61 

109 

113 

118 

20 

59 

64 

69 

120 

125 

130 

21 

66 

72 

77 

132 

137 

143 

22 

74 

80 

86 

144 

150 

156 

23 

82 

89 

95 

157 

163 

170 

24 

91 

98 

104 

171 

177 

184 

25 

100 

108 

114 

185 

191 

199 

26 

109 

118 

125 

199 

206 

215 

27 

119 

128 

136 

214 

222 

231 

28 

130 

139 

147 

230 

238 

247 

29 

140 

150 

159 

246 

255 

265 

30 

152 

162 

171 

263 

272 

282 

31 

163 

174 

184 

280 

290 

301 

32 

176 

187 

197 

298 

308 

319 

33 

188 

200 

210 

317 

327 

339 

34 

201 

214 

225 

335 

346 

359 

35 

215 

228 

239 

355 

366 

379 

36 

229 

243 

254 

375 

386 

400 

37 

243 

258 

270 

395 

407 

422 

38 

258 

273 

286 

416 

429 

444 

39 

274 

289 

302 

438 

451 

466 

40 

290 

305 

319 

460 

474 

489 

41 

306 

322 

336 

483 

497 

513 

42 

323 

340 

354 

506 

520 

537 

43 

340 

357 

372 

530 

545 

562 

44 

357 

376 

391 

554 

569 

588 

45 

375 

394 

410 

579 

595 

614 

46 

394 

413 

430 

604 

621 

640 

47 

413 

433 

450 

630 

647 

667 

48 

432 

453 

471 

656 

674 

695 

49 

452 

474 

492 

683 

701 

723 

50 

473 

495 

514 

710 

729 

751 
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TABLE II (Continued) 



0.99 

0.975 

0.95 

0.05 

0.025 

0.01 

51 

494 

516 

536 

738 

758 

780 

52 

515 

538 

558 

767 

787 

810 

53 

537 

561 

581 

796 

816 

840 

54 

559 

584 

605 

825 

846 

871 

55 

582 

607 

629 

855 

877 

902 

56 

605 

631 

653 

886 

908 

934 

57 

628 

655 

678 

917 

940 

967 

58 

652 

680 

703 

949 

972 

1000 

59 

677 

705 

729 

981 

1005 

1033 

60 

702 

731 

756 

1013 

1038 

1067 

61 

727 

757 

782 

1047 

1072 

1102 

62 

753 

784 

810 

1080 

1106 

1137 

63 

780 

811 

837 

1115 

1141 

1172 

64 

806 

838 

866 

1149 

1177 

1209 

65 

834 

866 

894 

1185 

1213 

1245 

66 

861 

895 

923 

1221 

1249 

1283 

67 

890 

924 

953 

1257 

1286 

1320 

68 

918 

953 

983 

1294 

1324 

1359 

69 

948 

983 

1014 

1331 

1362 

1397 

70 

977 

1014 

1045 

1369 

1400 

1437 

71 

1007 

1045 

1076 

1408 

1439 

1477 

72 

1038 

1076 

1108 

1447 

1479 

1517 

73 

1069 

1108 

1141 

1486 

1519 

1558 

74 

1100 

1140 

1174 

1526 

1560 

1600 

75 

1132 

1173 

1207 

1567 

1601 

1642 

76 

1165 

1206 

1241 

1608 

1643 

1684 

77 

1198 

1240 

1275 

1650 

1685 

1727 

78 

1231 

1274 

1310 

1692 

1728 

1771 

79 

1265 

1309 

1346 

1734 

1771 

1815 

80 

1299 

1344 

1382 

1777 

1815 

1860 

81 

1334 

1379 

1418 

1821 

1860 

1905 

82 

1369 

1415 

1455 

1865 

1905 

1951 

83 

1405 

1452 

1492 

1910 

1950 

1997 

84 

1441 

1489 

1530 

1955 

1996 

2044 

85 

1478 

1526 

1568 

2001 

2043 

2091 

86 

1515 

1564 

1606 

2048 

2090 

2139 

87 

1552 

1603 

1646 

2094 

2137 

2188 

88 

1590 

1642 

1685 

2142 

2185 

2237 

89 

1629 

1681 

1725 

2190 

2234 

2286 

90 

1668 

1721 

1766 

2238 

2283 

2336 
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TABLE II (Continued) 


0 



<x 








0.05 



91 

1707 

1761 

1807 

2287 

2333 

2387 

92 

1747 

1802 

1849 

2336 

2383 

2438 

93 

1787 

1843 

1891 

2386 

2434 

2490 

94 

1828 

1885 

1933 

2437 

2485 

2542 

95 

1870 

1927 

1976 

2488 

2537 

2594 

96 

1911 

1970 

2020 

2539 

2589 

2648 

97 

1954 

2013 

2064 

2591 

2642 

2701 

98 

1996 

2057 

2108 

2644 

2695 

2756 

99 

2040 

2101 

2153 

2697 

2749 

2810 

100 

2083 

2145 

2198 

2751 

2804 

2866 


♦Values of such that Prob[tj^ > tj^.^] = a. 


TABLE III. PERCENTILES OF THE CHI-SQUARE DISTRIBUTION.* 



Degrees 22 22 

of x,99 X.975 X.95 X .90 

Freedom 


1 

.00016 

, 

2 

.020 

. 

3 

.115 

. 

4 

.297 

. 

5 

.554 

. 

6 

.872 

1 . 

7 

1.24 

1 . 

8 

1.65 

2 . 

9 

2.09 

2 . 

10 

2,56 

3 , 

12 

3.57 

4 . 

14 

4.66 

5 . 

16 

5.81 

6 . 

18 

7.01 

8 . 

20 

8.26 

9 . 

22 

9.54 

11 . 

24 

10.9 

12 . 

26 

12.2 

13 . 

28 

13.6 

15 . 

30 

15.0 

16 . 

40 

22.1 

24 . 


00098 

,0039 

.0158 

051 

.103 

.211 

216 

.352 

.584 

484 

.711 

1.06 

831 

1.15 

1.61 

24 

1.64 

2.20 

69 

2.17 

2.83 

18 

2.73 

3.49 

70 

3.33 

4.17 

25 

3.94 

4.87 

40 

5.23 

6.30 

63 

6.57 

7.79 

91 

7.96 

9.31 

23 

9.39 

10.9 

59 

10.9 

12.4 

0 

12,3 

14.0 

4 

13,8 

15.7 

8 

15.4 

17.3 

3 

16.9 

18.9 

8 

18.5 

20,6 

4 

26.5 

29.0 


2 

X.IO 


2 

X.05 


2 

X.025 


2 

X.Ol 


2.71 

3.84 

5.02 

6.63 

4.61 

5.99 

7.38 

9.21 

6.25 

7.81 

9.35 

11.3 

7.78 

9.49 

11.1 

13.3 

9.24 

n.i 

12.8 

15.1 

10.6 

12.6 

14.4 

16.8 

12.0 

14.1 

16.0 

18.5 

13.4 

15.5 

17.5 

20.1 

14.7 

16.9 

19.0 

21.7 

16.0 

18.3 

20.5 

23.2 

18.5 

21.0 

23.3 

26.2 

21.1 

23.7 

26.1 

29.1 

23.5 

26.3 

28.8 

32.0 

26.0 

28.9 

31.5 

34.8 

28.4 

31.4 

34.2 

37.6 

30.8 

33.9 

36.8 

40.3 

33.2 

36.4 

39.4 

43.0 

35.6 

38.9 

41.9 

45.6 

37.9 

41.3 

44.5 

48.3 

40.3 

43.8 

47.0 

50.9 

51.8 

55.8 

59.3 

63.7 



TABLE III (Continued) 


Degrees 

of 

Freedom 


2 

X.99 


2 

X.975 


2 

X.95 


2 

X.90 


2 

X.IO 


2 

^.05 


2 

X.025 


2 

X.oi 
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TABLE IV. RESULTS OF RUN 10 MEASUREMENTS: TSI #1192 (TRIPLE SPLIT 
FILM). 


Block Means of the Velocity Components in 
Block Number Sensor Oriented Coordinate System 



v(A) fps 

v(B) fps 

v(C) fps 

1 

7.159 

0.315 

19.118 

2 

6.872 

1.232 

15.274 

3 

6.664 

12.136 

10.567 

4 

7.306 

11.930 

12.626 

5 

6.235 

5.927 

14.096 

6 

2.369 

3.524 

10.346 

7 

4.893 

1.206 

15.555 

8 

6.932 

-0.821 

23.212 

9 

6.692 

-1.322 

26.001 

10 

4.141 

2.949 

16.273 

n 

5.224 

1.336 

16.416 

12 

5.282 

3.082 

15.915 

13 

10.639 

14.491 

17.574 

14 

8,049 

3.766 

21.162 

15 

8.247 

6.195 

17.395 

16 

6.072 

0.441 

22.161 

17 

7,211 

0.205 

22.961 

18 

7.766 

0.980 

18.673 

19 

5.015 

-0.723 

18.898 

20 

5,915 

-0.143 

19.546 

21 

8,586 

2.636 

20.457 

22 

7.390 

3.662 

15.965 

23 

4.374 

6.836 

13.083 

24 

4,337 

-0.866 

19.375 

25 

3.472 

-1.763 

21.035 

26 

4.058 

-0.796 

20.246 

27 

6.941 

1.359 

18.494 

28 

7.503 

1.036 

20.262 

29 

9.147 

9.244 

15.837 

30 

7.301 

10.648 

11.629 

31 

8.605 

7.725 

12.869 

32 

7.905 

4.548 

16.498 

33 

4.665 

1.491 

16.278 

34 

6.505 

2.555 

16.601 

35 

7.209 

6.773 

12.281 

36 

5.946 

5.954 

11.613 

37 

9.551 

13.367 

13.501 

38 

7.881 

3.950 

14.587 
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TABLE IV (Continued) 


Block Number 

Block Means of the Velocity Components in 
the Sensor Oriented Coordinate System 

v(A) fps 

v(B) fps 

v(C) fps 

39 

5.730 

1.151 

16.624 

40 

1.979 

1.008 

11.641 

41 

2.102 

2.806 

11.127 

42 

0.468 

3.191 

8.227 

43 

1.599 

3.204 

7.521 

44 

7.733 

7.867 

12.748 

Sample Mean 

6.129 fps 

3.734 fps 

16.188 fps 

Number of Reverse 




Arrangements of 

504.0 

424.0 

584.0 

the Block Means 




Means, variances and covariances of the velocity components and 

temperature in the 

mean wind coordinate system: 



U 

V 

w 

T 

Z 

uv 

uw 

lifi. 

yw 

Yfi. 

s 


17.429 fps 
0.0 fps 
3.129 fps 
39.110 


15.778 (fps)^ 
-2.474 (fps)2 
-3.903 (fps)2 
-0.205 fps-^F 
18,029 (fps)2 
1.012 (fps)2 
0.513 fps-®F 
10.4 (fps)2 


-0.105 fps-®F 
0.455 (°F)2 
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TABLE V. RESULTS OF RUN 11 MEASUREMENTS: TSI #1192 (TRIPLE SPLIT 

FILM) - GILL ANEMOMETERS. 


Block 

Block Means 

of the Velocity Components 1n the 
Oriented Coordinate System 

Probe 

Number 

6111 Anemometers TSI #1192 
v'(x*) v'(x*) 

Gill Anemometers 
V(y*) 

TSI #1192 
v'(y*) 

1 

33.839 

31.774 

7.490 

8.561 

2 

27.660 

25.584 

7.645 

8.140 

3 

29.385 

26.814 

10.668 

11.404 

4 

22.272 

18.846 

18.131 

14.929 

5 

27.244 

23.009 

16.073 

16.952 

6 

26.837 

22.699 

16.747 

17.116 

7 

26.221 

22.869 

13.638 

14.326 

8 

20.513 

18.016 

12.177 

11.590 

9 

16.328 

14.732 

12.509 

9.598 

10 

24.719 

21.331 

17.402 

15.866 

n 

25.447 

22.354 

11.876 

12.642 

12 

16,918 

15.709 

9.262 

8.018 

13 

24.315 

20.883 

13.089 

13.879 

14 

26.210 

22.893 

12.069 

12.099 

15 

30.838 

27.496 

11.497 

12.394 

16 

32.487 

28,336 

14.951 

15.874 

17 

28.082 

25.038 

10.749 

11.963 

18 

26.172 

24.006 

8.004 

8.302 

19 

26.893 

25.376 

5.158 

6.399 

20 

24.310 

22.529 

7.533 

8.126 

21 

39.024 

37.139 

5.112 

6.206 

22 

42.911 

40.890 

4.986 

7.015 

23 

38.896 

36.741 

6.889 

8.387 

24 

40.385 

38.757 

2.037 

3.641 

25 

39.766 

38.061 

3.621 

5.018 

26 

37.026 

36.001 

3.353 

4,697 

27 

39.881 

38.289 

-4.187 

-3.475 

28 

30.765 

29.668 

-2.398 

-1.054 

29 

28.728 

27.510 

-1.806 

-1,240 

'30 

27.699 

26.396 

-3.908 

-3.883 

31 

27.268 

25.651 

1.018 

1.411 

32 

33.251 

31.501 

6.514 

7,734 

33 

28.653 

27.330 

2.609 

3.478 

34 

29.494 

28.187 

5.964 

7.159 

35 

31.008 

29.634 

5.831 

7.159 

36 

28.059 

26.900 

3.132 

4.153 

37 

24.640 

23.980 

-1.387 

-1.692 
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TABLE V (Continued) 


Block 

Block Means 

of the Velocity Components in the 
Oriented Coordinate System 

Probe 

Number Gill Anemometers TSI #1192 
v'(x*) v'(x*) 

Gill Anemometers 
v'{y*) 

TSI #1192 
v'{y*) 

38 

26.047 

25.087 

0.799 

1.076 

39 

25.754 

24.587 

3.797 

4.825 

40 

20.645 

19.657 

4.591 

4.639 

41 

22.335 

21.677 

-0.353 

-0.153 

42 

15.381 

15.310 

-1.259 

-1.701 

43 

17.851 

16.315 

-4.003 

-4.018 

44 

Number of 
Reverse 
Arrange- 
ments of 
the Block 

15.865 

15.173 

-2.678 

-2.659 

Means in 

Sensor 

Oriented 

Coordinate 

System 

512.0 

427.0 

766.0 

671.0 



Mean Wind Components in the Probe Oriented 

Coordinate System 


TSI #1192 

Gill Anemometers 

U' 

(fps) 25.93 

27.91 

V 

(fps) 6.70 

6.25 

T 

(°F) 33.03 



(deg.) -14,49 

-12.62 
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TABLE V (Continued) 


Means, Variances and Covariances of the Velocity Components 
in the Mean-Wind Coordinate System 



TSI #1192 

Gill Anemometers 

U (fps) 

26.778 

28.600 

V (fps) 

0.0 

0.0 

W (fps) 

5.920 


^ (fps)^ 

35.914 

2.464 

26.367 

0.319 

UW (fps)'^ 

-4.391 


ue fps-°F 

-0.743 


(fps)^ 

32.931 

22.680 

yw (fps)2 

0.050 


v0^ (fps-°F) 

-0.242 


wi (fps)^ 
we (fps-®F) 

25.316 

0.040 


^ (°F) 

0.367 



TABLE VI. RESULTS OF RUN 12 MEASUREMENTS: TSI #122 (DUAL SPLIT FILM) - GILL ANEMOMETERS. 


Block Number 

Block Means of the 

Velocity Components In the Probe 

Oriented Coordinate System 

Gill Anemometers 
v(x*) 

TSI #122 
v(x*) 

Gill Anemometers 
v(y*) 

TSI #122 
v(y*) 

TSI #122 
v(z*) 

1 

24.132 

23,535 

-3.656 

-0.149 

7.976 

2 

23,528 

23.560 

-1.027 

2.866 

8.096 

3 

28,992 

26.246 

-9.957 

-8.377 

6.793 

4 

32.051 

29.115 

-3.910 

-1.331 

7.950 

5 

32.846 

34.792 

-5.209 

-2.437 

7.064 

6 

28.414 

23,290 

-18.330 

-18.530 

7.804 

7 

30.508 

23.050 

-24.876 

-24.837 

7.244 

8 

26.462 

18.599 

-24.103 

-23.335 

5.787 

9 

22.803 

20,223 

-22.716 

-21.351 

6.405 

10 

22.515 

20,223 

-7.472 

-6.077 

6.498 

n 

32,998 

32.323 

2.572 

5.328 

8.175 

12 

35,420 

34.819 

1.880 

5.312 

8.931 

13 

19,690 

18.519 

-7.062 

-5.419 

5.366 

14 

22.059 

17.190 

-14.803 

-14.636 

5.289 

15 

16,780 

14.627 

-9.382 

-7.963 

4.142 

16 

16.681 

17.027 

-2.512 

0.433 

3.971 

17 

22.781 

21,644 

-0.704 

1 .207 

6.798 

18 

25.599 

24.852 

-0.144 

3.554 

8.288 

19 

36.042 

32.616 

-9.502 

-7.801 

6.364 

20 

23.357 

21.106 

-6.997 

-4.233 

7.121 

21 

18.160 

15.151 

-10.525 

-7.953 

4.729 

22 

23.236 

14,392 

-23.278 

-22.111 

6.049 

23 

19.137 

16.359 

-10.046 

-9.157 

6.888 

24 

17.285 

16.352 

-4.860 

-2.812 

4.520 

25 

23,981 

23.431 

-0.543 

2.057 

4.589 



TABLE VI (Continued) 



Block Means of the Velocity Components in the Probe 

Oriented Coordinate System 

Block Number 

Gill Anemometers 
v(x*) 

TSI #122 
v(x*) 

Gin Anemometers 
v(y*) 

TSI #122 
v(y*) 

TSI #122 
V(2*) 

26 

21.584 

20.232 

-7.411 

-5.219 

5.268 

27 

17.827 

11.484 

-15.543 

-14.812 

4.806 

28 

22.796 

15.873 

-20.170 

-19.084 

5.593 

29 

17.114 

15.384 

-5.864 

-3.682 

5.430 

30 

32.880 

30.078 

-7.292 

-5.822 

5.836 

31 

30.316 

28.692 

-4.863 

-1.108 

5.683 

32 

25.317 

23.307 

-5.360 

-2.556 

5.693 

33 

21.278 

17.926 

-10.252 

-8.780 

4.413 

34 

24.692 

19.871 

-15.139 

-14.429 

4.395 

35 

19.140 

16.372 

-9.616 

-7.810 

5.104 

36 

21.599 

18.726 

-10.101 

-7.680 

4.837 

37 

25,766 

21.717 

-10.178 

-6.826 

5.813 

38 

26.900 

25.464 

-3.067 

0.537 

5.895 

39 

39,215 

37.735 

0.576 

5.472 

7.868 

40 

29.653 

25.839 

-11.967 

-10.007 

6.523 

41 

24,603 

24.219 

-1.600 

1.987 

6.503 

42 

23.190 

19,188 

-13.438 

-11.190 

6.923 

43 

23.058 

17.221 

-16.795 

-14.658 

3.632 

44 

17,602 

14.790 

-11.366 

-9.221 

5.639 

Number of 

Reverse Arrange- 
ments of the 
Block Means 

524.0 

523.0 

514.0 

494.0 

601.0 



TABLE VI (Continued) 

Mean Velocity Components in the Probe Oriented Coordinate System 



^TSI #122 

c 

Gill Anemometers 

TSI #122 

Gill Anemometers 

U’ (fps) 

21.750 

24.772 

21.750 

24.772 

V (fps) 

-6.878 

-9.014 

-6.878 

-9.014 

W (fps) 

6.107 


6.107 


T (°F) 

41.609 


41.609 


6 (deg.) 

17.550 

19.995 

17.550 

19.995 


Means, Variances and Covariances of the Velocity Components in the Mean-Wind 

Coordinate System 



^TSI #122 

^Gill Anemometers 

TSI #122 

Gill Anemometers 

li (fps) 

22.812 

26.361 

22.812 

26.361 

V fps) 

0.0 

0.0 

0.0 

0.0 

W (fps) 

6.107 


6.107 


u2 (fps)2 

34.863 

27.760 

35.559 

28.320 

uv (fps); 

-0.216 

-3.393 

-0.116 

-3.239 

UW (fps)2 

-0.282 


-0.285 


iJ6 fps-°F 

-1.182 


-1.182 


— 

vw (fps)^ 

55.290 

35.451 

57.377 

36.160 

9.241 


9.412 


V 0 fps-°F 

-0.338 


-0.298 


(fps)^ 

20.873 


20.983 




TABLE VI (Continued) 



^TSI #122 

^Gill Anemometers 

TSI #122 

Gill Anemometers 

we fps-°F 

0.219 


0.227 


(°F)2 

1.452 


1.468 



^Block means removed. 




TABLE VII. RESULTS OF RUN 14 MEASUREMENTS: TSI #122 — GILL ANEMOMETERS. 


Block Number 

Block Means of the Velocity Components in the Probe 

Oriented Coordinate System 

Gill Anemometers 

V(>!*) 

TSI #122 
v(x*) 

Gill Anemometers 
v(y*) 

TSI #122 
v(y*) 

TSI #122 
v(z*) 

1 

27.628 

26.230 

14.520 

15.538 

6.729 

2 

27.143 

25.422 

14.621 

15.258 

6.440 

3 

26.295 

24.621 

13.140 

13.171 

5.488 

4 

18.198 

18.103 

6.020 

6.836 

4.866 

5 

25.149 

23.531 

10.516 

11.415 

6.182 

6 

30.664 

29.643 

14.508 

15.209 

7.989 

7 

25.101 

24.066 

10.748 

10.844 

5.784 

8 

26.115 

24.250 

16.188 

16.088 

6.887 

9 

31.156 

29.846 

10.589 

11.733 

6.228 

10 

26.985 

24,673 

-1.872 

-3.296 

7.295 

n 

31.578 

29.201 

18.478 

18.559 

7.719 

12 

35.638 

33.469 

19.658 

20.491 

8.862 

13 

28.971 

27.138 

9.979 

10.512 

7.365 

14 

26.242 

24.870 

12.074 

12.676 

6.943 

15 

23.638 

22.676 

10.260 

10.612 

5.101 

16 

23.057 

21.433 

13.697 

13.658 

5.845 

17 

31.742 

28.310 

20.198 

19.890 

7.806 

18 

33.644 

31.865 

17.692 

17.668 

8.193 

19 

26.085 

23.935 

15.668 

16.013 

6.009 

20 

25.228 

24.271 

10.891 

11.562 

6.729 

21 

25.775 

25.187 

2.402 

2.172 

6,951 

22 

24.197 

23.540 

7.788 

7.601 

7.842 

23 

28.916 

27.966 

12.146 

13.293 

7.141 

24 

23.690 

21.835 

-2.679 

-4.636 

6.452 

25 

24.901 

23.608 

-4.543 

-5.495 

7.405 

26 

20.647 

19.645 

-5.765 

-6.584 

7.069 
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TABLE VII (Continued) 


Block Number 


Block Means of the Velocity C omponents 1n the Probe Oriented Coordinate System 

Gill Anemometers TSf im Gill Anemomete"rs TSI #122 TSl #12: 

v()(*) v(x*) v(y*) v(y*) v(z*) 


27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

Number of 
Reverse Arrange- 
ments of the 
Block Mean 


25.743 

24.216 

-5.777 

-7.872 

19.740 

17.967 

-7.290 

-9.133 

26.619 

24.414 

12.687 

14.010 

26.126 

25.368 

4.391 

4.101 

27.358 

25.690 

-5.627 

-7.820 

22.217 

20.911 

-3.553 

-4.903 

18.328 

16.226 

-6.596 

-8.732 

27.728 

26.292 

2.808 

1.939 

28.495 

26,844 

10.775 

10.788 

28,370 

26.843 

8.828 

9.566 

25.417 

23.856 

11.189 

10.961 

34,486 

32.053 

19.427 

20.322 

45.617 

41.982 

25.594 

25.250 

28.522 

27,665 

13.154 

14.080 

37.873 

35.921 

16.237 

17.067 

32.370 

30.956 

4.584 

3.508 

21.294 

20.853 

4.196 

3.971 

19.579 

18.862 

2.896 

2.408 

472.0 

467.0 

558.0 

564.0 


7.205 

7.400 

7.710 

6.598 

7.291 

6.867 

8.258 

6.390 

6.362 

5.420 

7.837 

8.298 

9.975 

9.371 

9.696 

7.454 

3.768 

6.338 


375.0 



TABLE VII (Continued) 

Mean Velocity Components In the Probe Oriented Coordinate System 


— — 

StSI #122 

^Gill Anemometers 

TSI #122 

Gill Anemometers 

U' (fps) 

25.597 

27.142 

25.597 

27.142 

V* (fps) 

8.416 

8.519 

8.416 

8.519 

W (fps) 

7.035 


7.035 


T (“F) 

55.000 


55.000 


6 (deg.) 

-18.200 

-17.425 

-18.200 

-17.425 


Means , 

Variances and Covariances of the Velocity Components in the Mean-Wind 

Coordinate System 


§TSI #122 

§Gill Anemometers 

TSI #122 

Gill Anemometers 

U (fps) 

26.945 

28.448 

26.945 

28.448 

V (fps) 

0.0 

0,0 

0.0 

0.0 

W (fps) 

7.035 


7.035 


(fps)^ 

30.694 

27.325 

31.365 

28.088 

liv (fps)2 

4.043 

2.825 

4.231 

2.937 

uw (fps) 

-0.486 


-0.380 


U 0 (fps-°F) 

0.0 


0.0 


(fps)^ 

50.041 

35.886 

50.761 

36.391 

vw (fps)^ 

2.765 


2.776 


v0 (fps-°F) 

0.0 


0.0 



CT» 


TABLE VII (Continued) 



®TSI #122 

®Gill Anemometers 

TSI m2 

Gill Anemometers 

(fps)^ 

17.705 


17.758 


we (fps-°F) 

0.0 


0.0 


(°f2) 

0.0 


0.0 



"Block means removed. 



APPENDIX A 

THE RESPONSE CHARACTERISTICS OF THE GILL ANEMOMETERS 
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Sensor Response 

The Gill anemometer utilizes propellers, the shafts of which are 
parallel to each of the two perpendicular coordinate axes in the 
horizontal plane. The polystyrene propellers have four blades and 
are constructed to turn one revolution for each 31.7 cm of air 
passing the propeller. Each anemometer contains a miniature 
tachometer generator which is turned by the propeller and produces a vol- 
tage that is related to the respective wind components. Since the 
set of Gill anemometers were rotated with the TSI anemometers, they 
were oriented in such a fashion that one anemometer was on the 
average in the direction of the mean wind and the second anemometer 
perpendicular to the first one both in a horizontal plane. In this 
orientation both Gill anemometers were calibrated in the wind tunnel 
for velocities ranging from 10 fps to 70 fps and for various angles of 
attack varying from -40^ to +40^^ with respect to the mean-velocity 
direction. 

It is therefore suggested that this set of propeller anemometers 
should not be used when the angle of attack is larger than +40° or 
less than -40°. The velocities should not go beyond the range as 
specified. 

The velocity components can be obtained from the following 
empirical relations which fitted the calibration data very well for 
the above ranges in velocity and angle of attack: 
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and 


I 



0,069 


Ex 

(Cos B)®'^ 



= 0.085 Ey 


(A-1) 


(A-2) 


where Ex and Ey are the anemometer output voltages and 3 is the angle 
of attack. Since the angle B is not known at the outset, a simple 
iteration procedure is used to calculate the velocity component 
Vx- No single empirical relation was attempted to be obtained since the 
two propellers were always rotated in such a way so that one was 
along the mean wind and the other one normal to it. 

From the empirical equations obtained from the wind-tunnel calibra- 
tion data one can easily see that these anemometers deviate sys- 
tematically from the ideal cosine law and no reliance can be placed 
on just a single calibration with zero angle of attack. 

In turbulent conditions, there can be large differences between 
the indicated speed and the actual velocity. The differences arise 
from sensor sensitivity to relative wind direction, sensor dynamic 
characteristics, and to the data averaging procedure. The difference 
between measurements with various common types of anemometers may 
exceed 30 percent [62], 

The response time of the rotating type of anemometers such as cup 
anemometers and propeller type anemometers 1s always obtained by 
letting the anemometers accelerate from rest to some equilibrium 
speed of rotation depending on the wind velocity. The response of 
the anemometers which are allowed to decelerate from some speed of 
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rotation to rest is very difficult to realize in the wind tunnel. 

As a consequence the response time is based on the acceleration of 
the anemometers and the response for deceleration is assumed to be 
similar. However, the cup anemometers and propeller anemometers are 
known to accelerate faster than to decelerate with increasing and 
decreasing wind velocities respectively. As a result, the calculated 
mean velocity is estimated too high with respect to the actual velocity 
when the anemometers are used in turbulent winds. In addition, the 
output voltage of the anemometers shows a rectified undulation superimposed 
on the mean voltage. The frequency and magnitude of this "ripple" depends 
on the rate of rotation of the anemometer. If the propeller anemometers 
are used for the measurement of fluctuating velocities the digitized 
data will be affected by the presence of this undulation. 

The calibration data have indicated that the output of the propeller 
anemometer is related directly to the component of the velocity parallel 
to its shaft. Because of varying response characteristics when the 
propeller operates in a accelerating flow as compared to a decelerating 
flow, the mean velocity is overestimated (See results in Tables V, VI, 
and VII). Because of the limited response characteristics at frequencies 
higher than 1 hertz, the propeller anemometers underestimate the variances 
and covariances of the different velocity comopnents. This fact is very 
well illustrated if the spectra of the velocity components measured by 
the Gill anemometers are compared with those measured by the fast-response 
TSI anemometer (compare figure 18 with 21 and figure 19 with 22). 



APPENDIX B 

LISTING OF FORTRAN PROGRAMS 

I . DATPl 

II. TREND 

III. DATP2 

IV. TRANSFORM 
V, POWER SPECTRUM 
VI. CROSS SPECTRUM 
VII. SIMULATION 
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C DATA PROCESSING SYSTEM OF A SET OF FOUR DIGITIZED TIME SERIES 
C FOR CALCULATING ME AN, CORR ELAT I ON AND TRANSFORMATION MATRIX 
C U< If n=X-COMPONENTS VELOCITY OF SENSOR X»V(A) 

C U( 2fI)=:Y-CCMPON£NTS VELOCITY OF SFNSOP. 1,VIB) 

C U(3f n=Z-COMPONENTS VELOCITY GF SENSOR i,VtC) 

C U(4fl l=TeMPERATUR£ OF SENSOR 1 

C — — — 

DIMENSION U (4 1 8 192) fUBBK(44f 4) ,6ET AH44) f COVf 44f 4f 4) f DEG 1 44) 
DIMENSION UBKT< 44 ) i UBT0TL(4) ,TOTftEV( 8) ,TURl 3) , STD VI T( 44 , 8 t8 ) t 
lREVARR(44f8)fBKVAR<44f4f4) 


CALL DATA FROM MAGNETIC TAPE 
NOTBK=DATBK=NUMBER GF DATA IN EACH BLOCK 
NQBLK=IB=NUMB£R OF BLOCKS 
IX=NDTBK=NUHBER OF DATA IN EACH BLOCK 
IB=THE BLOCK NUMBER 
NTS-NUMBER OF TIME SERIES 


NTS=4 

NDTBK=8192 

NQBLK=44 

BLKNO=NGBLK 

IB=l 


nooooonn ooononnooonrjooosn 



o o r> o n 


IX=NOTBK 

5000 00 7 I^ltNDTBK 

7 PEADaCtlG) U(4,U ,U(2»n»U(3,I) 

10 F0RMAT(4A4) 


CALCULATE THE MEAN VALUE IN EACH BLOCK WHICH WE MIGHT CALL 
PIECEWISE MEAN 

UBBK( IB*J)-MEAN VALUE IN EACH BLOCK f 16= It 2t * . .t N08LK, J = l*2f.»4 


DO 11 J=1*NTS 

11 UBSKCIB, J)=0.0 
1 = 1 

12 DO 13 J==1,NTS 

13 U8BKnB,J) = UBBK(IB,J)+U<Jtn 
IF( I .EO. NDTBK) GO TO 14 
I»I+1 

GO TO 12 

14 0AT6K=NDTBK 
DO 15 J=1,NTS 

15 UBBK( IB«J) = UBBKUBf JI/0AT6K 

UBKTI IBI=SQRT(U86K( IB, U«*2*U8BKI I8t 2 » **2+U8BK U B, 3)**2) 

C CALCULATE THE TRANSFORMATION ANGLE BETAl UN EACH BLOCK 
C CALCULATE TNE MEAN SQUARES AND CORRELATIONS IN EACH BLOCK 
C COV=MEAN SQUARES AND CORRELATIONS IN A BLOCK 
C eKVAR=8L0CK VARIANCES AND COVARIANCES 

0ETAHIB) = ATAN(-1.22475<‘(UBBKCIB,2»-UBBK(IBt3l)/(UBBK(IBrl) 
l+UBBKt IB»2)+UBBKl IB,3)) I 
OEGII B)^{ 18 0.0/3. 14159)*BETAinBI 
DO 16 K=lfNTS 
DG 16 J=1 *NTS 


nonoon oc^nno 



o n o o o 


16 CUV( I8fK,J)=0,0 
1=1 

17 00 18 K=1,NTS 
DO 18 L=1»NTS 

18 CQV t IO»K,U = COV( IB,K,L M-U(K, n*U(Lf n 
IF(i -EQ. NOTBK) GO TO 191 

1 = 1+1 
GO TO 17 

191 DO 19 K=lfNT$ 

OQ 19 L=1,NTS 

19 COVUB,K,L)=COV{IB,K,L)/DATBK 
00 192 N=1*NTS 

BKV ARnB^NfNMC0Vn8,N,N)-UBBK<IB»N)*UBBK(IB,NJ 

192 ST0Vl TnB,N,N)=S0RT(A6S(BKVA«(IB,NtN> ) > 

IF (IB •EQ, N08LK) GO TO 3000 
IB=IB+l 

GO TO 5COO 
3000 CONTINUE 

KRITE 16,100) 

100 FORMAT( IHl, 1X,*6K N0»,2X,»MEAN VI A) • ,2X , • ME AN VI B ) * ,2X, • MEAN VIC)' 
1,3X, 'MAGNITUDE', 2X, 'MEAN TEMP • , 2X, • BLK ANG',3X,'ST0 OIV VA',2X, 
2'STD DIV V8',2X,'ST0 DIV VC',2X,'ST0 DIV T'//) 

DO 101 IB=1,N08LK 

101 WRITE 16,102) I8,UB6K( IB, n,UBBKlIB,2) ,U69K( IB, 3 ) , U8KT 11 B ) , U6 BK ( I P. , 
14),DEG( IB) , STDVITII B,1 ,1) ,STOV IT| IB,2 , 2 ) , STOV I T ( I B , 3 , 3 ) , STDVITUS, 
24,4) 

102 FORMAT (3X, I2,2X , F9 *3, 3 1 3X, F9 .3 ) , 2X, F9,3, 3X, F6 *2, 4( 3X, F9 . 3 ) ) 


HERE BEGINS THE CALCULATION OF THE MEAN OR WE MIGHT CALL IT THE 
GRAND OR ACCUMULATED MEAN, OR SAMPLE MEAN 
UBTOTU J) =ACCUMULATfn ME AN ,U=1 , 2,3 , 4 


n o o n o 



n r> Ci Ci or>oo 


DO 20 J=itNTS 

20 UBTOTLt J)=0,0 
J = 1 

21 DO 22 IB=ltKGPLK 

22 UBTOTL(J)=UBTGTL( J)+UBBKU6, J) 

IF(J .EQ, NTS) GO TO 24 

J=J + 1 
GO TO 21 

24 DO 23 J=1,NTS 

23 U8TOTU J)=U8T0TL<J)/BLKN0 
VT0TLi=SQ8T(UPT0TU 1 ) =^^2 + UBTOTL( 2 )*^2+UBT0T t ( 3 )*^2 ) 

WRITE (6, 104) UUBT0TLU)»J = 1»NTS),VT0TL1) 

104 FORMAT ( ///2X , • SA Mf>LE ME AN=» *4Fl2-6, 5Xf ’ TOT AL VEL. MAG .= • » F 12 . 6/ ) 


HERE BEGINS THE CALCULATIONS OF SAMPLE BETAS 
BETOLl=THE SAMPLE BETA 


BETOLl = ATAN(-l,22475*(U8TOTLl2l-U6TOTU3n/(UBTOTL( D+UBTQTLI2) 
l+UBTOTLOni 

6£OEG=(180-0/3.14159 J^BETQLI 
WRITE (6, 103) BETOLl fBEOEG 
10 3 FORMAT (/2Xf 'BETOL 1= • , Fia.8* 5Xt •BEDEG=SF 10, 51 


HERE CALCULATES THE NUMBER GF REVERSE ARRANGEMENTS FOR STU. DEV. 
IN EACH BLOCK 


NSTD=8 

NTEST=NC8LK-1 
DO 181 IB=1,N0BLK 
DO 181 r=i»3 
K = I+5 

181 STD VI T I IB,KfKj = STOVIT{I8,Itn 


or>no onon 



DO 182 IB=liN08LK 
STOVIT( ia,4,4)=UBKT( IB) 

STDVIK IB,5,5)=0£G( IB| 

DO 182 I=lf3 

182 5T0VI T(I8tI •! ) = U8BK( IBtU 
OQ 194 IB = l,aNOBLK 
DO 194 I=l,NSTO 

194 REVARfUIB,n = C.O 
18=1 

K=ia+l 

195 00 197 11*1, NSTD 
DO 197 JJ=K,NOBLK 

IF (STDVITt IB,U, in .GT. STDVIT ( JJ , II . 1 1) ) GO TO 196 

GO TO 197 _ 

196 REVARR( ID, lI)=REVARft( IB,ri) + 1.0 

197 CONTINUE 

IF (IB *£0* NTEST) GO TO 198 
I8*IB+1 

GO TO 195 

198 CONTINUE 

DO 190 J*1,NST0 
TOTREV( J)*0.0 
00 190 IB*l, NTEST 

190 TOTREVt J)*TOTREV( J)+ReVARRUB,J) 

WRITE (6,201 ) (T0TREV(JJ,J=1,NSTD) 

201 FORMAT(//10X, 'TOTAL NO- OF REV. ARR-*= ' , 8F1Q . 1 ) 

STOP 

END 



noon o ooooo 


:fr;^:jjt3{c:$:i{f J ffgQ:<c i{c <c * :{c »* ^ ijc <t *: <iS!! ^ !{{ :<£:^£ :*s 


THIS PROGRAM IS US6D TG APPLY THE MOVING-AVERAGE AND DIFFERENC 
HIGH-PASS FILTER TC THE ORIGINAL TIME SERIES. 


DIMENSION U(4»X6384) ,UBBK( 4 ♦ 8 192 ) 


N08LK=NUMBER tlF BLOCKS CONSEQUENTLY FILTERED. 

NDATA=NUMBER OF DATA POINTS TO BE USED AT THE START OF FILTERING 
PROCESS. 


NDATA= 16384 
NOBLK=43 
NDTBK=NDATA/2 
NHALF*N0T8K/2 
NTS*4 
18*1 

DAT6K=:N0TfiK 
NEW=NOTBK+l 
DO 1 I=l,NOATA 

1 READ(10,2) U(4,I),U{l,n,U<2,I),U(3»U 

2 F0RMATI4A4) 

GO TO 2000 

1000 DO 3 I*N6W,N0ATA 

3 READ(io,2j u(4, n ,u{ un ,UI2,n,UI3,n 

2000 DO 4 J=1,NTS 

4 U8BK( J,1 1=0.0 
DO 5 J=1,NTS 
DO 5 I=1,NDT8K 

5 UBBK( Jil)=UBBKUf 1 )+U(J»Il 
DO 6 1=1, NTS 

6 UBBK( I, 1 l = UeBK{ I, 1 1/DAT3K 


OU OO uo uO 



DO 7 1=1, NTS 
J=NOTBK+l 

7 U8BK( I,U = UP8K( I , 1) + (U U , J )-UU , 1 )> /OATBK 
DD 8 J=l,NTS 

DO 8 I=2,NDTBK 
K = I-1 
L=NOTBK+I 

8 UBBK* Jt I)=UBBK( J,K)+(U(J,L|-U(JtI))/DATBK 
DO 9 K=1,NTS 

DO 9 I=1,NDT6K 

J=NHALF4-M 

9 U(K,n = U(K,J)~UB8K(K,T) 

DO 10 I=l,NDT8K 

10 wRiTE(ii,2» u(4,n,u(i,n*ui2fn*u(3,n 
DO 11 LL=1,NTS 

DO 11 I = 1,NDTBK 
M=I 

J=NDT BK+M 

11 U(LL,I)=UCLL,J) 

IF (18 .eQ* NOBLK) GO TO 5000 
I8=I84‘1 
GO TO 1000 
5000 CONTINUE 

WRITE(6,5C50) IB 

5050 FORMAT (20X, 'NUMBER OF BLOCKS COMPLETED*,! 
STOP 
END 
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DATA PROCESSING SYSTEM CALCULATING VARIANCES AND COVARIANCES C 

BASED ON THE DATA READ FROM THE TAPE WITH EQUALLY-WEIGHTED C 

RUNNING MEAN APPLIED. SAMPLE MEAN AND BETA ANGLE ARE READ FROM C 
PREVIOUS PROGRAM- THE PRINTED VALUES ARE THE STATISTICAL VALUES C 
IN THE MEAN WIND COORDINATE SYSTEM- C 

THE MEANING OF THE CODE NAMES MAY BE READ FROM OATPl- C 

DIMENSION U(4»8192) , COV ( 43 ,4,4) ,UBTOTL 1 4 ) * SC0V(4 ,4) ,TUB( 3) 
DIMENSION El(4,4)fUHEAN(4) ,COVMENI4f 4) tTBINTi4) 

DIMENSICN X8AR(4) *XSUM(4) 

NTS=4 

NOT 8K=8192 

NOBLK=43 

6LKN0=N0BLK 

DAT8K=NCTBK 

IB=1 

UBTOTLI 1)= 6-128789 
UBTOTL(2)= 3-733873 
UBTQTU3I=16. 187836 
UBTOTH 4 1 = 39.109543 
6ETCL1= 0.529701 


SINCE THE MEAN HAS BEEN REMOVED SO WE CALCULATE THE MEAN SQUARE 


C 

C 


o o 



non noon ono 


VALUF WHICH SHOULO BE THE VAIUE CLOSE TO THE VARIANCES AND C 

CiJVARTANCES 



5000 00 7 I=l,NDTRK 

7 READtlOaO) U(4,n,U(Un»U(2#I),U<3,l) 

10 F0RMAT(AA4> 



HERE THE SMALL NEAR ZERO MEANS ARE CALCULATED AND SUBTRACTED FROM 
THE DATAS IN THE BLOCK. C 

DO 81 1=1, NTS 

81 XSUM(I»=0.0 
DO 82 J=1,NTS 
DO 82 K=l,NOT6K 

82 XSUM(J)=XSUM(J) + U(J,K) 

00 83 l=l,NTS 

83 XBAR( I)=XSUH{ IJ/NOTBK 
DO 84 JaltNTS 
DO 84 I=1,NDTBK 

84 U(Jf I) = U( Jf I)-XBAR( J) 


BLOCK VARIANCES AND COVARIANCES ARE CALCULATED IN THE FOLLOWING 


DO 16 K=1,NTS 
DO lb J=1,NTS 

16 C0V(1B,K, J)=0.0 
I=i 

17 DO 18 K=1 ,NTS 
DO 18 L=1,NTS 

18 COV ( IB,K,U = C0V( I6,K,L)+UCK,n«UlL,l) 
IF! I .EO. NDTBK) GO TO 191 

1 = 1+1 


n n o 



nonoo nr?nor> 


GO TO 17 

191 DO 19 K=1,NTS 
DO 19 L=1,NTS 

19 COV( IB,K,L)=CaVUB,K,U/DAT8K 
IF (16 .EQ, NOBLK3 GO TO 3000 
ia=I 8+1 
GO TO 5000 
3000 CONTINUE 


here begins the calculation of the sample variances and covariances 
FOR the whole length OF THE DATA 
SCOV=SAMPLE VARIANCES AND COVARIANCES 


DO 25 K^lfNTS 
DO 25 J^l.NTS 

25 SCOV(K* Jl=0-0 

26 DO 27 K=lfNTS 
DO 27 LslfNTS 

27 SCOV(K,U-SCOV(K,U+COV(IB,K,U 
IF(I8 .EO. ND8LK) GO TO 28 

I 8=IB+1 
GO TO 26 

28 DO 29 K=1,NTS 
DO 29 L=l,NTS 

29 SCDV(K»L)=SCQV(K,U/BLKNO 


HERE BEGINS THE TRANSFORMATION OF MATRIX 

El(-=i,4J IS the matrix TRANSFORMATION FROM THE SENSOR DIRECTION 
TO THE MEAN WIND DIRECTION 


DO 9? 1=1, NTS 


or>ooo ono 
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Uf'iEAN ( 1 ) = 0*0 
DO 43 J-l,MTS 
EH I,J)=0.0 
43 COVMENn, J)=0,0 

EU 1, U=0.5773 5*CGStBeT0Ll) 

El( 1, 2)=0,57735*SIN(BET0tl) 

El(l,3)=-0. 81650 

El(2. 1)=0*57735*COS<BET 0Ll)-0 .7071 1 *S I N( B ETOLl ) 

EH 2*2)=‘0-57735#SIN(BETOLI)+0.70711*COS{6ETOL1) 

£1(2,31=0.40824 

El( 3, 11=0. 57735*C0S(8ET0LU+0.7071l#SIK!( BET OLD 
E1(3,2)=0.57735#SIN(BETOL l)-0.70711*caS( BETOLl) 

E1(3,3)=C.4C824 
£1(4,41 = 1.0 

UMEAN(I) IS THE MEAN WIND VELOCITY COMPONENTS & TEMPERATURE C 

COVMEN( I, J ns THE VARIANCES AND COVARIANCES OF WIND TURBULENCE IN 
THE MEAN WIND DIRECTION. C 

Q 

1=1 

46 DO 45 J*1,NTS 

UMEANd )=UMEAN( I) +E 1 ( J, D*UBTOTU J) 

00 45 K=l,NTS 
DO 45 L=1 ,NTS 

45 COVMENf I,J)=COVMEN( I , J ) •*■£ U K , D*E1 ( L , J DSCOV( K , L ) 

IF (I .EQ. NTS ) GG TO 47 

1 = D1 

GO TO 46 

47 CONTINUE 


WRITTEN STATEMENTS 


u u c 



i^RITE(6,105) 

105 FOR^^ATUHlfAX, 'BEFORE TRANSFQRMAT ION ■ t27X MEAN WIND DIRECTIONS’/) 
WRITei6,106> 

106 F0RMATI2X, 'CORRELATION NO* •» 3X, • SAMPLE VAR. £ COVAR, • ,10Xt • VAR.E 
ICOVAR.' ,5X# 'VELOCITIES'/) 

DO 107 K=lrNTS 
DO 107 L=1,NTS 

107 ViRITE(6,108) K, L, SCQV (K» L) ,COVMEN CK,L ),UM EANIK ) 

108 FORMAT ( 7X , 12 ,3X , 12 , 8X,Fl 0.3 ,15X»F1 0.3 ,8X ,F8 *3 ) 
VT0TLl^SQRTiUBT0TUl)«*2+L©T0TU2)#«2+UBT0TH3)##2) 

DO 32 I-lt3 

32 TBINTd )=SQRT(SCOV( I, I) )/VTOTLl 
WRITE (6,207) C TB INT (I ) , 1=1 ,3 ) 

207 F0RMAT(///2X, 'TURBULENT I NTE NSITIES=' ,3F10.6 ) 

DO 113 1=1,3 

113 TU8(n=SQRT(C0VMEN( I, I) )/VTOTLl 
WRITE(6,114) (TUB(1 >,1=1,3) 

114 FORMAT (///2X, 'TURBULENT INTENSITY IN THE MEAN WIND 0|RECTI0N=', 
13F8.4) 

STOP 

END 
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— transform the data from the sensor oriented direction to the 

MEAN WIND direction FOR SPECTRUM CALCULATIONS-1192 £ 1193 

THE MEANING OF EACH CODE NAME MAY BE SEEN FROM DATPl. 


DIMENSION U(4f8192) , V(4 i 8192 ) t El ( 4t4 ) 

NTS=4 
NDT6K=8192 
NOCHK^l 
NOBLK=38 
RETOLl* 0,6165 

TRANSFORMATIONS C 

DO 3 I=lfNTS 

00 3 J»ltNTS g 

3 eu I,J)*0,0 

EKlf U=0*57735«COS(BETOLU 
EU li2)=0,57735*SIN(B6T0Lll 
EU I, 3J=-C,81650 

El(2»l)=0.57735*CCSieETntl)-0.70711*SIN(BETOLn 
El( 2, 2)^0. 57735*5 IN<BETOLU+0.70711*COS(BETOLI) 

EU 2.3)=0.4C824 

EU 3t 1)=0. 57735* COS (BET OL I U-0,7071 1 *SIN( BETOLU 
EU 3,2)=0-57735*SIN(BETOLl)-O.7O711*CnS(BETGLl) 

EU3,3) = 0,40824 
EU4, 41=1.0 
1000 DO 1 I=ltNDTBK 

1 readuOtIo) u( 4, n,u( If 11 tU(2,n,u(3f n 
10 FQRMAT(4A41 

DO 2 I=lfNTS 
DO 2 J-l I AOT3K 

2 V( I ,J)=0.0 


u u o u u 



1=1 

46 DO 45 K=ltNOTBK 
DO 45 4=liNTS 

45 V( I ,K)=V( !,K) + EUJ»n^U(JtK> 

IF { I •FQ. NTSI GO TO 47 
1 = 1+1 
GO TO 46 

47 CONTINUE 

DO 9 I=l,NOTBK 

9 WRITElll.il) V(4tI),V{l,I).V<2,I),V<3.n 
11 FORMAT (4A4I 

IF (NOCHK .EQ. NOBLK) GO TO 2000 
N0CHK=N0CHK+1 
GO TO 1000 
2000 CONTINUE 

WRITEI6.999) NOCHK 

999 F0RMAT(5X. 'NUMBER OF BLOCKS TRANSFORMED' , 15 ) 
STOP 
ENO 
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POWER $P EC TR UM*##******#** *«***^*^^*f'*’''*^ 


THIS PROGRAM! CALCULATES THE POWER SPECTRAL DENSITY ESTIMATES 
INCLUDING TEE DATA TAPERING, FAST EGURIER TRANSFORM AND SMOOTH! 
TECHNICUES. 


CIMENSICN U (8 192 ), SP 11 4097), SP2( 40971, PS11(9A),SPP( 4097) 
DIMENSION FLECM94) ,FNE0(94) ,FMeO<94) 

DIMENSION F(8192 ) ,0( 820) ,H( 820) 

DIMENSICN C (2049) ,X (8192), Z(2048),V( 8192) 

DIMENSION FRE0(94) ,FLOG(94),PNORM(94) 

COMPLEXES E,X,2 


NDTBK=NUMB£R OF DATA POINTS IN EACH BLOCKED TIME SERIES 
CAT BK^NCTBK 

N0BLK*^U«8eP OF BLOCKS USED IN THE SPECTRAL ESTIMATE. 

MAVE=NUM8ER OF SPECTRAL VALUES AS A RESULT OF COMBINED SMOOTHING 
NTS=NUM6EP OF TIME SERIES TO BE PRECESSEO 
NnwiO=NUM6ER OF DIFFERENT DATA WINDOWS TO BE TAPERED 
NCK=CCNTRCL NUMBER USED FOR NUMBER OF DATA WINDOWS 
KOUNT=TFE CCNTROL NUMBER CEPENCS ON TEE TIME SERIES NEEDED. 

0 IVIS=MCNDC“1/10 NUMBER OF TAPERING POINTS 


N0TBK=8192 

DATBK=NDT8K 

NHALF=NCTBK/2 

HALFN=NHALF 

MAVE=94 

NP=94 

CN=^MAVF 

CC=CN 

C.U=CN 


ooooor>oor»oo nf^o 



Dr=0.CC5 

N08LK=43 

PLKNQ=NC8LK 

NTS=2 

KCK = 1 

NGWID=1 

KGUNT=1 

PI=3.1^15^ 

KGDNO=a20 
DIVIS=MGDKC 
NHF=NHALF+1 
ViRITE(6,20l ) 

201 F0RMAT(lHl,50Xf 'CCSINE TAPERING WINDOW*//) 

C c 

C D(I) AND H(I) ARE THE END TAPERING FUNCTION IN COSINE TAPER DATA C 

C WINDOW* F(I) IS THE NUMBER ONE DATA WINDOW USED FOR TAPERING C 

C C 

00 12 I*UNOTBK 
P*I 

12 F(I) = l.O-((P-(DAT0K-'l.O)/2.vO)/HDATEK + i.O)/2.D))**2 
DO 14 I=1,MC0N0 
P = I 

140(1 )=C.5*( 1.0-CCS(PI*P/DIV1S) ) 

DO 16 I=1,MCDN0 
P=I 

16 h(I .0-CCS (P (DATeK-P)/CIVIS ) ) 


f Q 

C THIS PART CALCULATES THE FRl DUENC I E S tDF IS THE LGwPST C 

C CNYC IS THE NYQU! ST FREQUENCY. C ANOW 0=8 A NO n I D TH C 

C FLhD=EREQ=FNEQ=FNEQ=DIFFERENT P EPRES tAT 'iT ICNS Op FPEOUENCY VALUES 

r C 

i>= 1.0/ ( CAreK'•^OT ) 
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QNYQ=i,C/(2.0=f^OT) 

6ANDWC=CN=«'DF 
C1=DF 
0 2=D1 
D3=D2 

CG 66 1=1,8 
P = I 

66 FLGaU)=P*OF 
DO 61 1=9,38 
P = I 

t=P-6,0 

61 FLEQ(I)= (4.0*(Q-i,0) + 2.5)*DF 
DO 62 1=39,62 

G=I 

P=Q-5.0 

62 FLEQm = (l6.0*(P-25.0) +9.5)«=C1 
DG 63 1=63, 86 

P = i 

R=P-6*0 

63 FLEQl n= (6A.0*(R-49.0 J +33. 5)^2 
no 64 1=81,94 

R=I 

SS=R-6.0 

64FLEOU)= (256.0*( SS-73.C) +129.5)^D3 
DO 65 1=1, NP 
FREQC IJ=FLEQ< I ) 

FNEOt I ) = FLEQ( I ) 

65 FM£Q{ J) = FLEQ( I ) 


DATA READING FROM TAPE 

KCUNT=TFE CCNTROL V4LUF5 USED TO PF/D CIFFEPFNT DATA SfTS 


u U <>; o 



na lOA 1=1, NHF 
104 SP2< I )=C.O 
4000 10=1 

00 TO ti ,2,3) , KOLiNT 

1 READ { 10,11 ) U 
11 FQR^ATU4) 

GO TO ICCC 

2 REWIND 1C 

4 RcAD(iO,2ll U 
21 FQKMAT(AX,A4) 

GO TO 1000 

3 REWIND 10 

5 R£AD( 1C, 31) U 
31 F0RMA7(12X,A4) 

1000 CONTINUE 

C C 

C FILTERING OF THE TIME SERIES MTH MOVING AVERAGE ALREADY APPLIED C 

C TO ADJUST TFc TIMF SERIES TC NEAR ZERO MEAN, SMALL MEAN VALUES C 

C ARE AGAIN CALCULATED AND SUBTRACTED FRCM THE VALUES, OF DATA POINTS 

C IN EACH BLOCK. 

C C 

XSUM=0.0 
DC 81 I = 1,NDTBK 
81 XSJM =um +XSUM 
XbAR=XSlM/OATBK 
nc 82 I = 1 tN 0TBK 
8? U I ) = Un ) -XBAR 


C C 

C STATEMENTS LSCD TL NULTIPLY DIFfFRFNT WINGCWS C 

C 

C:0 ru ( 14*1,148) ,NCWID 
144 OCJ 18 i=l,MCDNO 



lb t(n=u(n^cn) 
DO 17 I=l,MCDNO 
K=I+7372 

17 0(K)=t(Kl^H(I) 
GO TO 500 

145 OG 18 K=l»NDTeK 

18 U(K)=U*K)=S=FfK) 
530 CONTINUE 


C c 

C CALCULATE THE FOURIER CQ6FFIENTS BY USING FAST FOURIER TRANSFORM C 

c c 

DO 99 K=l,NOT6K 
99 V(K) = 0,0 

CALL FFT (UfV,XJ 

C c 

C CALCULATE THE SPECTRAL VALUES TOTALING 4097 POINTS. THESE VALUPS 

C ARE ACJLSTED OY A SCALE FACTOR DUE TC COSINE TAPERING. LOOP 103 C 
C IS TO SU'iMING StGMENTLY ALL THE BLOCKED SPECTRAL VALUES. LOOP 109 

C rS TO AVERAGE THEM. SPP( H ARE THE VALUES SAME AS SP2CI). C 

C C 

DO 102 1=1, NHF 


102 SPIU )=2.C>«'rT’»'(ReAL< X ( I ))**2+AINAG( X ( I ) )/DATBK 
IF (NC^IO ,NE. i) GC TO 93 

CO 101 K=l, NHF 

101 SPI(K) =(1.0/0,8751*SRl (KI 
98 DO 103 1 = 1, NHF 

103 SP2 ( 1 ) = SP2< I J+SPl ( I ) 

IF (IB ,£0. NC6LK J GC TO 6000 
I B= I B + 1 

GG TG (1 ,4,5) , KCUM 
6000 CONTINUE 

CC 109 T=l,NHr 



SP2 ( I ) = SP2 ( I)/3LKN0 
DO 121 1=1, KHF 
121 SPPd J = SP2(n 

WR I TP (6, 222) lSP2t I ),I = l,JMHF) 

222. FQRMAT(5X,iCF10*4 1 

C FOR FR6QUENCV SMOOTH ING ,TH£ ESTIMATED SPECTRUM MAY BE CONSIDERED 

C AS REPRESENTING THE MIDPCINT CF THE FREQUENCY INTERVAL 

C SINCE ONLY HALF NUMBER OF THE POINTS ARE UNIQUE AFTER TRASFORM C 

C THE SPECTRAL AVERAGE IS PERFORMED BASED UPON THE TOTAL OF 4096 C 

C VALUES, CCNSEQUENTLY 94 SMGGTHEO SPECTRAL VALUES ARE CALCULATED, 

DC 22 1=9,94 

22 PSIK I) = C,0 
DO 71 K=2,9 
J=K-i 

71 PS11< J) = SP2IK) 

L = 10 
MA=1B 
M4=4 
K = 9 

24 DC 23 I=L,MA 

23 PSil(K)=PSll(K)+SPPn ) 

IF (K ,EC- 38) GO TO 25 
K=K+1 

L=L+M4 
MA= MA+M4 
GO TO 24 

25 CONTIFUE 

PC 42 J=9,3S 
42 PSIK J)=PSIU J) /4.C 
L = 130 
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M16=16 
K = 39 

52 OG 53 I = 

53 PSlHK) = PSlHKi+SPPn) 
IF (K .EC* 62) GC TO 54 
K = K + i 

L=L+M16 
P8^HB+M16 
GU TQ 5 2 

54 CONTINUE 

DO 43 J=3<3f62 
43 PSU( J) = PS1HJ)/16,C 
1=514 
NC=577 
^:64 = 64 
K = 63 

55 DC 56 I=LiMC 

56 PSiK K) = PSll(K) + SPPn ) 
IF (K .FC. 86) GC TO 57 
K = KU 

L=L+M64 
NC=PC+M64 
GO TO 55 

57 CONTiNUr 

DO 44 J=63f86 
4^ PSUI J) = PSll(J)/64.G 
I. =2 05C 
f'': U= ? 'i 0 5 
^256=256 

5 b L‘0 5S I=Lt^o: 
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5s Psiiu)==Psu(K)+spp(n 

IF (K .EC, 94) GC TO 60 
K=K+1 
L=L+M256 
P.D=MD+M256 
GO TO 56 
60 CGNT INUE 

no 45 J=a?.94 
45 PSIK J)=PS11(J)/256.C 
DO 27 LL=1,NP 

27 FL0G( LL) = AtCG(FREC(LUJ 

DO 28 1 = 1, NP 

28 PNORMU »=FKEO(n«PSil(n 


ALL TEE WRITTEN STATEf^ENT 


IF (NOWID .EQ. 2 .AND. KCUNT .EC. 1) GO TO 32 
GO TO 36 

32 WRITE(6,33I 

33 FQRMATt IHlt SOX, 'WINDOW NUMBER 1'//) 

36 CONTINUE 

KRITEC6tlC6) 

106 FORMAT ( lOX ,* SMOOTHED SPECTRUM* ,5X, ' FREQUENCY • ,5 X ,• FREO . MULTPlD 
1SPECTRUN*,5X,'LCG OF FREQ.'/J 

DC lOT I=1,NP 

107 WRITE (6, ice) PSIK n,FR£Q( n,PNORM( n,f LGC< I J 

108 FORMAT 1 1 4X , FIO .4 , 8X , F8 .4, 14X,F 10. At 8X,F 8.4 ) 


THE CCNTROLLED VALUE GF NTS WILL OBTAIN THE CALCULATIONS OF 
DIFFERENT NUMBER GF TIME SERIES. THE CONTR(>LL^D VALUE OF NCK 
WILL CALCLLATE DIFFERENT DATA WINDCWS. 


noooo noo 
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IF (KCONT .EQ. MTS) GO TO 3C0C 
KCUM=:KCUNT + i 
on 110 1 = 1, NHF 

110 SP2 (1 ) = Q*0 
Gc in 4000 

3000 CONTINUE 

IF (NCWIO .EQ. NCK) GO TO 37 

NQriIO=NCWID+l 

no in I=1,NHF 

111 SP2(I)=0.0 
REWIND 1C 
K0UNT=1 

GC TO 4000 
37 CONTINUE 
STOP 
END 

SU380UTINE FFT (U,V,X) 

IMPLICIT REAL*4( A-H,P-Z1 

GIMHNSICN C(2049)*X<8192),Z(2048),lJ(81^2) ,V( 8192) 

COMPLEXES E,X,Z 

C 

C N = TCTAL NC. CF DATA 
C NS = STAGE NO • 

C NJ = TOTAL NO. CF ACCITIDN ( AND SUPTRACTION ) STEPS IN EACH STAGE 

C NU = TOTAL NO. OF ADDITIONS ( OR SUET. ) IN ONE AUDITION STEP 

C lA = STARTING NO. OF EACH ADDITION STEP 
C NT = ENDING NC. Of EACH ADDITTCN STEP 
C J = AUDITICN STEP NC. IN EACH STAGE 

C I , KA = SUBSCRIPT OF NEW X IN EACH STAGE 

C INVERT=1 IS USED TO OBTAIN THE INVERSE FOURIER TRANSFORM. FOR ANY 

CTHER NUMRERS OF INVERT INDICATES FCURIEP TRANSFORM. 



INVERT=3 
N=8192 
NN=N/ 2 
KM=NN/2 
KK=NN+1 
■ km=nm+i 

CALL COSINE (N,NN,N^,C) 
DC 12 r=l,N 

12 xn )=c^PLX(U( n ,v(i) ) 
NS =1 
NJ-NN 
NU=1 

I M=NJ/2 
IA=l 
NT=NU 

DO 71 L = l,Nf' 

LA=L+NM 
71 ZiU=X(LA) 


N A=MM 
NB = 1 
KA^KK 

DC 2 4 = 1 tNJ 

na=na-nl 

. NB=iMB + NU 

IF (J-NI) 41,42,43 

41 E=CMPLX(C(NA),-C(NB)) 

CO rn 42 

43 ir=CMPLX(-CtNA),-C(NB)) 

42 CO P I=IA,NT 

IF U.LC.M) GO TC 81 




IC=I-NH 
IB=I+NN 
K A^KA-fl 
K=KA^-^U 

X(KA)=Z(ICI+X(IBJ 
x(K )-(znc)'X{i8jHe 
GO TQ 9 

81 1C=MM-I 
le=lC+NK 
KA=KA-1 
K=KA-NU 

IF ( J.EC.M ) GO TC 82 
X(KAJ = (X( IC)-X(l6n*F 

GG TO 83 

82 'X(KA)aXnC)-XnB) 

83 X(K) = X( IC)4X( IB * 

9 CONTINUE 

I A=IA+NU 
NT=NT+NU 

IF (J.NE.M) GO TO 22 
NA=MM+M’ 

NB=1-NU 
KA=NN 
GO TQ 2 
22 KA=K 
2 CONTINUE 

C 

IF (NJ-EQ-.2) GO TO 11 

NU=2*«NS 

NS=NS+l 

NJ=NI 

GO TO i 



C 
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11 CO .:.L I = 1,NK 
I B= I+NN 
C = X( ! )4X( IB ) 

X{ £B)=X( I )-X( 18) 

21 xn )=E 

C 

THE NEXT EIGHT CARDS ARE USED CNLY IF INVERSE FOURIER TRANSFORM IS 
NECESSARY. IN THIS CASE, THIS SUBROUTINE IS USED TO OBTAIN A REAL 
time series um instead cf using for spectral calculations. 


IF (INVERT .EQ. 1) GO TO 27 
GO TO 28 

27 CG 23 1=^1, N 

23 xii )=xn )/N 

U m = real (X(D) 

CG 25 I»2 *N 

25 U (I) «REAL(X(N-I+2)) 

28 CONTINUE 
RETURN 
END 

SUBROUTINE COSINE (N,NN,NM,C) 
DIMENSiCN C (20A9) 

TN=NN 

ANG=3 .1A1S927/TN 
I NTP^N/8 
CS=CQS (ANG) 

SN=SIN <ANG) 

C(1 . 

I I==Nr' + l 
C( 1 i )=0. 
or 39 J^l, IMF 
J I = JH 


u u 



JJ=NM+2-J 

C( JI)=C(J)*CS-C(JJJ’!'SN 
JM=NM+ 1-J 

39 C(Jf<) = C( J)^SN+C(JJI*CS 

RETURN 

END 
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c C 

C THIS PRfGR^H CALCULATES THE PCWER SPECTRUM, CO-SPECTRUM, QUAORATU- 

C PH SPECTRUM, CROSS SPECTRUM, PHASE ANGLE AND CDHERENCY FUNCTION C 

C Cy GOING THROUGH DATA TAPERING, SPECTRAL SMOOTHING TECHNIQUES. C 

C c 

D IMENSIGN U(8192) ,V( 6192) , SPA(A097> ,SPB<4097) ,SPC (4097) ,SPD(4097), 


1SPU4Q9 7),SP2(4C9 7) , CSP ( 4097 ) ,OSP ( 4097 ), FREQ( 941 , FLOG ( 94) 
OIMENSICN AA(1) , BB ( I ), URE (4097) ,U £M 14097 ) ,VR E( 4097 ), V IM ( 4097) 
OIMENSIGN F(8192) ,0(820) ,H(820) 

DIMENSION C (2049) ,X( 8192) ,2(2048) 

DIMENS I CN PSll (94) ,PS22(94) ,CS12 (94 ) , QS 12 (94 ) 

DIMENSION AMP12(94)tCQHl2(94),PHS(94) ,DEGPS(94) 

0 IMENSICN P ACM! (94) ,PN0M2(94),CSN0M(94) ,0SN0M(94) 

OIMENSICN FL60(94) ,FNEQ(94) ,FMeO(94) 

D IMENSICN CRQS(94 ),CRSNM(94),CRS( 10) 

COMPLEX AA,BB 
COMPLEX*a E,X,Z 


C C 

C EXPLANATION OF THE DIFFERENT SYMBCLS USED IN THIS PROGRAM CAN BE C 

C CBTAINEC FPCM THE POWER SPECTRUM PROGRAM* C 

C c 


NDT8K=8i92 

DATBK=NCTPK 

NHALF = NCTeK/2 

HALFN=NHALF 

r,HF=NHALF + l 

.MAV£=94 

NP=94 

CN=MAVE 

NOBLK=43 

hlknc=ncblk 



yi 
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KT5 = 1 
f^CK=l 
NOW 10 = 1 
K CU NT = i 
DT=0,OOS 
PI=3, 14155 
^'CDN0=820 
CIVIS=MCDKC 
WRI TE(6»201) 

201 FGRMAT<1H1,50X, •C05iN£ TAPERING WINDOW*//) 


WINDOW CALCULATIONS 


DO 12 I=l,NOraK 

P=I 

12 F(I) = 1.0-MP-(CAT8K-l,0)/2.0)/Hr)ATBK+l.O)/2*0n’i'*2 - 

no 14 I= 1 ,MC 0 N 0 S 

P=I 

14 D( I) = C. 5 *( 1 . 0 -CCS (RI*P/CIVISn 
DO 16 I* 1 ,MCDN 0 
P=I 

16 HU J=C. 5 #( 1 , 0 -CUS (PI <=(DATBK-PJ/OI VIS) ) 

THIS PART CALCULATES THE FREQUENCIES 

C 

DF=l,C/( CATBK*DT ) 

CNYy=l,C/(2,0^DT) 

3ANDWC=CN*DF 
C 1 = 0 F 
C2 = D1 
03=02 

DO 66 I = 1 f 8 


non 
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P = I 

66 flBQi I ) = P*DF 
Ca 61 

PPP=! 

0=PPP-6,0 

61 FLfcO(I>=: (4*0#JO-1.0) + 

DO 62 1-39,62 

SS=I 

P=SS-6,0 

62 Fi.E0Cn = (16-0*(P-25*0) +9,5)*01 
DC 63 1=63, 86 

nLP = I 
B=RLP-6,C 

63 FLEGm= (64.0*(R-49.0) >33.5)*D2 
DO 64 1=87,94 

PQR=I 

SS=PGR-6.0 

64 FL£0(I)* (256.0*t SS-73.0J +129*5)#03 
DO 65 1=1, NP 

FR£Q( n=FL£C( 1 I 
rN£Q( I )=FL£C(1) 

65 FM£g( I)=FLECi I ) 

DO 27 LL=l,NP 

27 FLGG(LL) = ALOG ( F PE G ( LL ) ) 


DATA READING FROM TAPE 


DO 10^ 1=1, NHF 
SPi f I ) = G.O 

sp 2 (n=o.o 
c 3P ( n=c,o 

104 t)SP( I )=C ,0 


O O 



ri o n 


4010 I3~l 

GO TO ( 1,2,3 J I KCUM 

C CROSS SPECTFUM FOR U AND W C 

1 DO 10 I=l,NCTBK 

I 0 READ (10,11) Un ) ,V (1) 

II FQRMAT(AX,A4,4X,AA) 

GC TO 1000 

C CROSS SPECTRUM FOR U AND V C 

2 REi^INC 10 

4 GC 20 I^1,NCTBK 

20 READ* 10,21 ) U( I) , V( I ) 

21 fURMAT(4X,2A4) 

GG TO 1000 

C CROSS SPECTRUM FCP T AND H C 

3 REW IND 10 

5 DO 30 1=1,NDT6K 

30 PFAD( 10 , 31 ) un ) , vn ) 

31 FCRMAT(A4,8X,A4) 

1000 CONTINUF 


NEAR ZERO MEAN VALUES ARE REMOVED FROM THE FILTERED TIME SERIES. 


XSOMU^C. 0 
XSUMV=0.0 
Cl si T=l,NnTFK 
XSUMU =U(I) +XSUMU 
XSUMV ^V( I) +XSUMV 

xbaru=xsumu/oathk 

XPAHV=XSUMV/DATDK 
OlJ B2 I = 1,NCTBK 
U( I ) =U( ! ) -XDARU 
V U ) = V ( I) -XBARV 


o n o 



O O o o o o o 


APPLY niFFFRENT DATA WINDOWS TC THE DATA POINTS. 


GO TO (1A^,145),NGWID 

144 CO 15 I = 1#I^CCNC 
V( I )-Y(I )*0( I) 

15 Un ) = U( I )*D( I ) 

DO 17 I = l,^^CDNO 

K=I+7372 

V(K) = V(KJ*HU ) 

17 U(K)=U{K)*H(I) 

GO TO 5CC 

145 DO 18 K=1,NCTBK 
V(K»=V(K)«F(K) 

13 U(K)*U(K)*F(K) 

500 CONTINUE 


THE FAST FOURIER TRANSFORM IS CALCULATFO 8Y STORING UU) IN THE 
REAL PART ANO Vd) IN THE IMAGINARY PART. 


CALL FFT (U,V,X) 

AA(1)=CRPIX( REAL(XUI)»0.0) 

BBt 1)=CPPLX(AIMAGIX(I) J ,0.0) 

UBE ( I ) = HEAL( A At 1 ) ) 

UIHl I )==AIRAG( AA( I J ) 

VRtIl) = REAL(BBI 1) ) 

VIMIl )=AIMA0(BB1 1 ) ) 

DC 444 1=2, NHF 
K=NOTPK“I -♦-2 

UREil) = REALi( y ( I ) + X (K* )/2.0) 

UIMd ) =AIMAGH X ( I ) - X (K) )/2-C) 

VREU) =AIRAGM X ( I ) + X IK) )/2.0) 


o n o o o o 



on noon 


444 ViM(I) =-PEALM X ( I ) - X (KJ )/ 2 * 0 ) 

DO 102 1=1, NHF 

SPAt I )=2.0*CT*(URE(I)«=»2 ♦ U IM ( U*#2 ) /DATBK 
Span )=2.0^DT^(VRE( 1)^*2 + VI>'n)^<‘2)/DATBK 
SPC( U=2 .0*CT*(UPE( U=»VRe(l) + U I«( I )«VI MU ) ) /DATBK 

102 SPUU)=2.0^DT*UjRe( I)=i‘VIM(n - UIMi n^VREi n )/OAT8K 
IF (NQWID -NE. 1) GG TG 93 

DC 101 X=l,RHf 

SPA{K) = « 1.0/0- 875) ^^'SP A ( K) 

SP3 (K)==ll.0/C.375)*SPBIK ) 

SPC(K) =(1.0/0.875 )*SPC(K) 

101 SPO(K)=(1,0/0.875)*SPD( K> 

98 DO 103 1=1, NHF 

SPIU ) = SPl( I) + SPA(I) 

SP 2 ( 1 3=sp2( n + spen ) 

CSP( I )=CSP( I) SPC( I ) 

103 QSPU) = QSP(n + SPOU) 

IE U8 .EQ. NOBLK) 60 TO 6000 
18=18+1 

GO TO (1,4,5) , KOUNT 
6000 CONTINUE 

FOR SEGMENT SMOOTHING, SUMMING NUMBER OF THE CORRESPONDING POINTS 
IN EACH DESIGNATED BLOCKS. ALSO AVERAGED THE SUMMED VALUES. 


DO 109 1=1, NHF 
SPl U ) = SP1( n/BLKNC 
SP2tI )=SP2(I)/BLKNG 
CSP( I )=CSP( 1 )/BLKNO 
109 QSP ( I) = GSP( n/8LKNC 


ECk frfcuency smoothing, the estimated spectrum may be cdnsidfred 


o r-) n n 



ri o o o o 


AS REPRESENTING THE HiDPGINT Cf THE FREQUENCY INTERVAL 
SINCE GNLY HALF THUMPER OF THE POINTS ARE UNIQUE AFTER 
THE SPECTRAL AVERAGE IS PERFCFMEC BASEC UPON THE TOTAL 
POINTS AND TOTAL CF POINTS ARE OfiTAIEO 


DU 113 1=1.10 

113 CRSU) = SQRTtCSPd + QSP(n=**2l 

DO aaee 1 = 1.10 

8 388 WRITE(6,S999 J SP 1 ( 1 1 , S P2 ( IJ . CRS ( I ) . C SP( I» , QSP (I ) 
9999 FORMATdOX .5F16.41 
DO 22 1=9, NF 
PS111I)=0.0 
PS22U)*C,0 
CS12( n=c.o 

22 GS12I 1 ) = 0,0 
OQ 71 K=2,9 
J=K-1 

PSll ( J)=SPl (K) 

PS22U) = SP2(K) 

CS12( J)=CSP(K ) 

n QS12( 

L = 10 
RA=13 
H4=4 
K = 9 

24 CO 23 i = L,MA 

PS11(*<I=PS11 (KI +SFI ( I ) 

PS22< KI = PS22( K) +5P2(T ) 

CSL2(K) = CS12(K)+CSPn I 

23 QS12(K)=0S1?(K)+QSPU) 

IF IK .EQ, 38) GO TO 25 
K=K + 1 


TRASFCRM 
OF 4096 


o o n o 



L^L+M4 
MA=MA+M4 
GO TO 24 
25 CONTIME 

DO 42 J=9,38 
PSIK J)=PSlUJ)/4.0 
PS22( J)=PS22 ( J)/4 .0 
CS12( J)=CS121J) /4-C 

42 GS12( J)=QS12( J )/4 ,0 
L=130 

MB=145 

f^l6=l6 

K=39 

52 DO 53 I=L,MB 

PSU(K) = PS1HK) + SP1( I) 
PS22(KI=PS22(K)+SP2(r) 
CS12{K)=CS12IK J+CSP( U 

53 OS12(KJ*OSl2(K)+OSPm 
IF (K .EQ. 62) GO TO 54 
K=K+i 

L=L+M16 

+M16 

GO TO 52 

54 CONTINUE 

DO 43 J=39,62 
PS11(J) = PS11U)/16,0 
PS22{ J)=PS22( J)/X6.C 
CS12(J )=CS12( J)/16.0 

43 QS12( J)=QS12( J)/16.0 
L = 514 

MC=577 

?^ 64=04 


ro 

o 


PO 



K^63 

35 00 56 I = L,^^C 

PS11{K)=PS11<K)+SP1(I I 
PS221K >=PS22(K) + SP2( i I 
CS12(K)=CS12(K)+CSP< I > 

56 0S12<KJ=QS12(K)+CSP(I) 
IF IK .EG. 86) GO TC 57 
K=K + 1 

L=L+M6A 
wC=MC+H6A 
GO TO 55 

57 CONTINUE 

00 44 J=t3,66 
PSIK J) = P$U{J)/64.C 
PS22( J)»PS22( J)/64.0 
CS12( J)*CS12( JJ/64.C 

44 GS12( J)=CS12( J )/64.C 
L=205C 
MO* 2305 
M256=256 
K=87 

5ft CO 59 I=L ,MD 

PSi 1 {K) = PS11(K )+SPl( I ) 
f‘S22(K)*PS22(K) + SP2(l) 
CSi2(K)=CS12(K)+CSPtI) 

5P QSi2( K)=QS12(K)+QSPM ) 
IF (K .EQ. S4) GO TO 60 
K = K + 1 
L*L+M256 
MU=‘-'0 + iV25fc 
GO JO 58 

6 0 CGNTINUF 



ro 

o 


CO 



n o o 


Du ^5 J=e7,S4 
P5U( J)=PSiliJ)/256.G 
PS22(JMPS22(JJ/256.0 
CSJ.2( JJ=CS12(JJ/256.C 
45 CSi2U) = GSi2t J)/256.0 
00 114 1=1, NP 

114 CRUSC n = SGRT(CS12n + QS12II>^^2) 

DO 265 1 = 1, NP 

AMP12(I)=CS12(n=*'CS12(I) + 0S12(I)«'QS12(I) 

265 PHS(n = + ATAN{QS 12 (Z)/CSl 2 Cn 3 
00 266 1=1, NP 

COH12(n=AMP12{I) / (PSll(U*PS22n )) 

266 DEGPS< n=18C.0=f'PHS(n/(2.0^PI) 

DO 28 1=1, NP 

PNQMUl)=FREQtI )=^PS1UI) 

PNQM2 (I ) = FR£Qn)=«'PS22U ) 

CRSnM(I) =FRE0{I)^'CRCS(I) 

CSNOMt I )=FREQ( 1)^0512(1) 

28 QSNONn) = FREQn)#0S12n) 


ALL THE WRITTEN STATEMENT 


IF {.NOWlD .EO. 2 -AND. KOUNT ,EQ. 1) GO TG 32 
GO TO 36 

32 WRITE (6, 33) 

33 FQRMATdHl, SOX, 'WINDOW NUMSER 1'//) 

36 CONTINUE 

I'RITE: (6,106) 

106 F0R^AT(5X, • SMCD PCWEP SPFCT 1*,5X,»SML0 POWER SPECT 2 ' , 5X , • COSREC * 
1 , 5X , * GUAOSPEC * ,5X COHERENCY* ,5X, 'PHASE! DEG) *, 5X , 'EKE 0* , 2X, * CROSS 
2SPEC* /) 
on 107 1 = 1, NP 


u u o 



o o o o 


107 re (6,103) psii (i>,PS 22 n ) , csi 2 ( u,csizn ),cohi 2 (i),degps(I), 

IFREQ( n,CP0S( I ) ^ 

108 fORMAT(10X,F10,4,10X,F10*4f6X,Fe*4,6X»FB*4,6X,F3.4,6X,F8*2,4X,F«.4 

i,jX,F8*3) 

HRl Tt(6,2C3) 

203 FORMAT (//5X, * FREO PJPD SPEC i**5X,*FREC MUTLO SPEC 2S5X,*FREQ MUT 
ILD COSPEC* ,5Xt* FREO RUTLD OUACSPEC ' ,3X , * LOG OF FREQS 3X, *FRFO MULT 
2CRQSS* / ) 

DO 204 1=1, NP ^ , 

204 WRl T6 (6,205 ) PNOMl ( I ) ,PNGM2C I),CSNOMC I ), QSNOH( I ) , FLOG! I ) ,CRSNM( I ) 

2Q5 FORMAT( ICX ,F 10 . 4, 15X, FI 0*4, 10X,F1 0-4, lOX ,F1 0- 4, lOX , F8 .4 , lOX ,F8-4) 

NTS CCNTRCLS THE NUMBER OF TIME SERIES BEING CALCULATED. NCK 
CONTROLS THE NUMBER OF DATA WINDOWS TO BE USED. 

IF (KCUNT -EQ. NTS) GO TG 3000 S 

KCUNT = KCU NT + 1 
DO lie 1=1, NHF 
SPl (I )=0.0 
SP2(l )=0,0 
eSPd J = C,0 
no QSP( i ) = C.O 
GO TO 400C 
3000 CONTINUE 

IF (NCWID -EG- NCK) GO TO 37 
NGWID=NDWID+1 
DO 111 1=1, NHF 
SPl I I )=0,0 
SP2U ) = 0.0 
CSP( I ) = 0-0 
111 OSPU )=C.O 
REWINO 10 


o n o o 



K0JN7=1 
GO TO 4CQ0 
37 CDNTIME 
STOP 
END 

SUBROUTINE FFT (UtV»X> 

IHPLiCIT PEAL«4 (A'HtP-Z) 

COMPLEXES E,X,Z 

DIMHNSICN C 120^+9) ,X( 81921, Z<2048),U( 8192 3 ,Vt 8192) 

C 

C N = TOTAL NG. OF DATA 
C NS = STAGE NC . 

C NJ = TOTAL NO, OF ADDITION ( AND SUBTRACTION ) STEPS IN EACH STAGE 

C hU = TOTAL NC. OF ADDITIONS t CR SUET. ) IN ONE ADDITION STEP ro 

C U = STARTING NO. OF EACH ADDITION STEP § 

C NT = ENDING NO. OF EACH AOOITICN STEP 

C J = AODITIGN STEP NO, IN EACH STAGE 

C I , KA = SUBSCRIPT OF NEW X IN EACH STAGE 

C C 

INVFRT^3 

N=8192 

AN=N/2 

NM=NN/2 

KK=NN+1 

DO 12 1=1, N 

12 XU j =CIWPL Xtut n , VU ) ) 

HM=NN+l 

CALL CCSINF (N,NN,NM,C) 

NS=l 
NJ = NN 


r> o n 



NU=l 

1 NI=NJ/2 
IA=1 

NT = NU 

DC 71 L=l,N^^ 
LA=L+NM 
71 7(Lj=X(LAJ 


KA = MM 
N8=l 
KA = KK 


DO 2 J=1,NJ 
NA=NA-NU 
NB=!^B + NU 

IF (J-N!) Al,42,43 

41 E=CMPLX(C(NA),-C(NB)» 
GO TO 42 

43 E = CMPLX(-C{NA),-UNB) } 

42 DO 9 I=iA,NT 

IF (J .LE.M) GO TC 81 

KA=KA+1 

K=KA+AU 

X(KA)=Z(IC)+xn6) 

X(K ) = (ZUC)-X( IB))#E 

GG TO 9 
81 IC=MP-l 
I8=IC+NN 
K A=KA-i 
K = KA-MJ 




n o r> 


If ( J.EC.M) GO TO 82 
X{KA) = {XUC)-X( iS)>*E 
GO TO 83 

82 X(K4)=X{iC)-X(lB) 

8? X{K) = X( IC) + X(IBJ 

« CONTINUE 
T A=IA + NU 
NT=NT+NU 

IF (J.NE.Nn GO TO 22 
N'A=MM+NU 
N9= 1-NU 
KA = NN 
GO TO 2 
22 KA=K 
2 CONTINUE 


IF (NJ. E0.2I GO TO 11 
NU=2^«NS 

NS=NS+1 
NJ=NI 
GO TO 1 

ll DG 21 1=1, NN 
I 6= i+NN 

f_=xm+x(iBi 

xii&)=x(i)-xnB) 

21 X( I )=f 


THE Ft?LLOWING FIGHT STATEMENTS ARE USED ONLY IF INVERSE FOURIER 
TRANSFORM IS NECESSARY IN THIS FFT PROGRAM. C 

C 

IF (INVERT ,FO. IJ GC^TC 27 
GC Tfl 2 8 


n n 



Z1 DO 23 1^1, N 
23 X (I ) = X( n/N 

U ( 1) ^ SEAL (Xtl) ) 

DO 25 !=2fN 

25 U (1) =REALlXtN-T+2)) 

28 co^riNue 
RETURN 

ENi) 

SUaROUTINE CQSIME *NtNN,Nf*,C) 
DIMENSION C(2049) 

TN=NN 

ANG=3.1415S27/TN 
I NT B=N/a 
CS=COS (AKG) 

SN^SIN lANG) 

Cfl ) = !• 

1 I-NM+1 
C( IIJ=0, 

DO 39 J=l, INTB 
J I = J+ 1 
JJ=NM+2-J 

C( J1 ) = C(J J*CS-C(JJ )»SN 
JM=NM+l-J 

3V C t JM)=C( J)*SN + CIJ J) 

^ RETURN 
END 



simulation^ '^'*=^****<'*=«^*»**^**^«***^**#*5><'C 


C c 

C THIS PROGRAM SIMULATES A TIME SERIES RY USING THE VON KARMAN SEMI- 

C EMPIRICAL SPECTRUM BY USING FAST FOURIER TRANSFORM METHOD. C 

C THE SAMPLING PATE IS 20 SAMPLES PER SECOND TO CREATE THE 32768 C 

C TOTAL number OF DATA POINTS FOR APPROXIMATELY OF 1638 SECONDS C 

C ^ 

DIMENSION 4(32768,1,1), M(3), INV(8192), $(8192), PP(32768) 

COMPLEX A 
UBAR = 40.0 

COEF = SQRT<32768. 0*20.0/2.0) 

C THE VON KARMAN SEMI*EMPIRICAL SPECTRUM C 

C BOTH GAUSS AND HARM ARE SCIENTIFIC SUBROUTINE IN THE ISM PACKAGE 

C C o 


DO 100 I = 2,16384 
AI = I 

FREQ =» (AI-1.0)/1638.4 
SPECT = (57.6*UBAR)/ 

1 ( ( I -0+70. 78* (< 36 0.0*FREQ/UBAR I **2I )**( 5. 0/6.0) ) 
H = SORT(SPECT) 

CALL GAUSS ( 19 , 0.707, 0.0, XI) 

CALL GAUSS (13, 0.707, 0.0, ETA) 

XREAL = C0EF*H*XI 

XIMAG = CCEF*H*£TA 

Ad, 1,1) = CMPLX(XREAL, XIMAG) 

100 CONTINUE 

DO 200 J = 1,16383 
K = 16385 + J 
L = 16385 - J 

A(K,1,1) = CGNJG( A(L,l,in 
200 CONTINUE 



A(ia*l) = CMPLX (32768. 0=^UBARf0.0> 
AU6385»ia) = CMPLX(0.0,Q.0l 
M(ll = 15 
M(2 ) = C 
y(3) = 0 

CALL HARM(A,M,lNVfSt-l»IFERR) 

DO *iOO 1=1,32768 
400 PP( I)=REAL(An,l,in 
DO 300 N = 1,32768 
300 WKiTE<10,50) PP(N) 

50 FORMAT (A4) 

WRITE (6,3000) ( PP (I ), 1=1 ,32768,32) 
3000 FORMAT(5X,12FIO*2I 
STOP 

END .. .. 
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